Method and system for determining magnetic susceptibility distribution

ABSTRACT

Systems and methods for determining a distribution map of susceptibility property of an object are provided. The method may include one or more of the following operations. A phase diagram corresponding to a magnetic resonance (MR) signal of the object may be obtained. A preliminary field map may be determined based on the phase diagram. Preliminary error limiting information associated with the preliminary field map may be obtained. A preliminary distribution map of susceptibility property of the object may be determined based on the preliminary field map and the preliminary error limiting information. An iteration process including at least one iteration may be performed to determine a target distribution map of susceptibility property of the object.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority to Chinese Patent Application No.201810323015.X, filed on Apr. 11, 2018, the contents of which are herebyincorporated by reference.

TECHNICAL FIELD

This disclosure generally relates to a method and system for magneticresonance imaging, and more particularly, to a method and system fordetermining a magnetic susceptibility distribution map of an object.

BACKGROUND

Almost all material can become magnetized in a magnetic field and thendisplay magnetic effects similar to magnets (e.g. generating a localmagnetic field around the material). Magnetic resonance imaging (MRI)uses externally imposed magnetic fields (including a main magnet and agradient magnet) and radio frequency waves to generate images of anobject. The object is magnetized in the externally imposed magneticfields and generates local magnetic fields. The local magnetic fieldssuperimpose with the externally imposed magnetic fields, and may affectthe signal of MRI.

Susceptibility property (or referred to as magnetic susceptibility ormagnetic susceptibility property) is a measure of the extent of amaterial becomes magnetized in a magnetic field, or in other words, theintensity of a local magnetic field generated by the material (ifmagnetized). The susceptibility property of an object can be used toderive information such as the content of iron, the degree ofcalcification, and/or the blood oxygen level in tissue or an organ ofthe object.

Quantitative Susceptibility Mapping (QSM) is an emerging technique inMRI for determining the susceptibility property of an object. A generalprocess of QSM technology includes obtaining one or more sets ofmagnetic resonance phase information or phase images, generating a localfield map based on the phase information, and determining a distributionmap of the magnetic susceptibility based on the local field map and aformulated relationship between the local field map and the distributionmap of the magnetic susceptibility. The local field map can beapproximated to a convolution of a dipole kernel and the distributionmap of the magnetic susceptibility. Accordingly, the distribution map ofthe magnetic susceptibility can be determined based on the local fieldmap and an inverse dipole kernel. Unlike the dipole kernel, the inversedipole kernel contains singularity values, thus making the inverseproblem mathematically ill-posed and leading to apparent streakingartifacts in the solution results.

The current QSM techniques are facing at least two problems. The firstproblem is that there may be computational related streaking artifactsin the QSM maps due to the ill-posed nature of the inverse dipolekernel. Superficially, the streaking artifacts occur mainly due to theexistence of the singularity values in the inverse dipole kernel. Thusthe streaking artifacts can be reduced when the degree of singularity isreduced (which is the basic principle of the truncated k-space QSMmethods). However, the truncated k-space QSM method(s) can change theinverse dipole kernel, and may lead to underestimated susceptibilityvalues. And even the artifact is reduced, the remnant streaks usuallyremains visible in the final distribution maps. The co-localization of apriori information and the susceptibility distribution may reduce thestreaking artifacts by adding regularization to the QSM. However, thediscrepancy between the a priori information and the actualsusceptibility distribution may lead to either over- orunder-regularization, affecting the accuracy of the results. Also fortissues or regions with strong susceptibility variance, such asair-tissue boundaries or large bleedings in tissues, relatively strongstreaking artifacts may be produced, affecting other regions ofinterest.

The second problem arises from anatomical a priori information that isinvolved in solving optimizing QSM as mentioned. The a prioriinformation may provide a reference for the spatial characteristics ofthe susceptibility distribution (e.g., spatial continuity within tissue,and/or boundaries between different tissues). With such spatialinformation, one can impose a certain restriction in the form of, e.g.,regularization terms (or Bayesian statistics), in solving the ill-posedinverse dipole problem. However, the obtained a priori information maybe inaccurate, and the inaccurate a priori information may generateconvergence errors in QSM. Usually, the a priori information isestimated by, e.g., the magnitude images of the same scan, or an initialestimate of the susceptibility distribution map. However, there can benon-negligible mismatch between the estimate of the a priori informationand the actual spatial continuity of the tissue and/or boundariesbetween different tissues. Also, the estimate of the a prioriinformation is considered constant during the calculation process ofcurrent QSM techniques, hence any mismatch in the a priori informationestimation may lead to convergence errors of the calculations.

Therefore, it is desirable to provide a method and system fordetermining a distribution of susceptibility property of the object withreduced convergence errors and suppressed streaking artifacts.

SUMMARY

According to an aspect of the present disclosure, a method fordetermining a distribution map of susceptibility property of an objectis provided. The method may be implemented on a computing device havingat least one storage device storing a set of instructions and at leastone processor in communication with the at least one storage device. Themethod may include obtaining a phase diagram corresponding to a magneticresonance (MR) signal of the object and determining a preliminary fieldmap based on the phase diagram. The method may further include obtainingpreliminary error limiting information associated with the preliminaryfield map and determining, based on the preliminary field map and thepreliminary error limiting information, a preliminary distribution mapof susceptibility property of the object. The method may includeperforming an iteration process including at least one iteration. Eachof the at least one iteration may include updating error limitinginformation in a current iteration based on a distribution map ofsusceptibility property of the object updated in a preceding iteration;updating the distribution map of susceptibility property of the objectin the current iteration based on the error limiting information updatedin the current iteration; and performing a comparison based on theupdated distribution map in the current iteration and the preliminaryfield map. The method may further include designating the updateddistribution map of susceptibility property of the object as a targetdistribution map of susceptibility property of the object.

In some embodiments, the performing a comparison based on the updateddistribution map in the current iteration and the preliminary field mapmay include generating an estimated field map based on the updateddistribution map of susceptibility property of the object in the currentiteration; comparing the estimated field map with the preliminary fieldmap; and determining whether the updated distribution map in the currentiteration satisfies a condition based on a result of the comparisonbetween the estimated field map and the preliminary field map.

In some embodiments, the comparing the estimated field map with thepreliminary field map may include determining a difference field mapbetween the estimated field map and the preliminary field map; anddetermining an evaluation value based on the difference field map. Thedetermining whether the updated distribution map in the currentiteration satisfies a condition may include comparing the evaluationvalue with a threshold.

In some embodiments, the error limiting information may include an apriori map, and the obtaining preliminary error limiting information mayinclude determining at least one of a magnitude diagram, a phasediagram, a susceptibility weighted imaging (SWI) map, or a R2* (T2*) mapbased on the MR signal; and determining a preliminary a priori map basedon the at least one of the magnitude diagram, the phase diagram, thesusceptibility weighted imaging (SWI) map, or the R2* (T2*) map.

In some embodiments, the error limiting information may include an apriori map, and the updating error limiting information in a currentiteration based on the distribution map of susceptibility property ofthe object updated in a preceding iteration may include determining afeature map in the current iteration based at least in part on thedistribution map of susceptibility property of the object updated in thepreceding iteration; and updating the a priori map based on the featuremap.

In some embodiments, the distribution map of susceptibility property ofthe object updated in the current iteration may include a plurality ofpixels, each pixel of the distribution map of susceptibility property ofthe object in the current iteration having a pixel value. The featuremap may include a plurality of pixels having a one-to-one correspondencewith the plurality of pixels of the distribution map of susceptibilityproperty of the object. The determining a feature map in the currentiteration based on the distribution map of susceptibility property ofthe object updated in the preceding iteration may include, for each ofthe plurality of pixels in the feature map, determining an evaluationvalue between a pixel value of a first pixel and a pixel value of asecond pixel adjacent to the first pixel in the distribution map ofsusceptibility property of the object in the preceding iteration, anddesignating the evaluation value as a pixel value of a pixel in thefeature map of the current iteration that corresponds to the first pixelin the distribution map of susceptibility property of the object.

In some embodiments, the determining an evaluation value between a pixelvalue of a first pixel and a pixel value of a second pixel next to thefirst pixel in the distribution map of susceptibility property of theobject in the preceding iteration may include determining the evaluationvalue between the pixel value of the first pixel and the pixel value ofthe second pixel next to the first pixel in a reference direction, thereference direction including at least one of an X direction, a Ydirection, or a Z direction.

In some embodiments, the error limiting information may include an apriori map. The updating error limiting information in a currentiteration based on the distribution map of susceptibility property ofthe object in a preceding iteration may include determining an a priorifield map in the current iteration based on the distribution map ofsusceptibility property of the object in the preceding iteration, andupdating the a priori map in the current iteration based on the a priorifield map in the current iteration.

In some embodiments, the determining an a priori field map in thecurrent iteration based on the distribution map of susceptibilityproperty of the object in the preceding iteration may includedetermining the a priori field map in the current iteration byperforming a forward dipole kernel function on the distribution map ofsusceptibility property of the object in the preceding iteration.

In some embodiments, the error limiting information may include anartifact estimation map. The updating error limiting information in acurrent iteration based on the distribution map of susceptibilityproperty of the object in a preceding iteration may include generatingan estimated field map in the current iteration based on the updateddistribution map of susceptibility property of the object in the currentiteration, determining a difference field map between the estimatedfield map in the current iteration and the preliminary field map, andupdating the artifact estimation map in the current iteration based onthe difference field map.

In some embodiments, the generating an estimated field map in thecurrent iteration based on the updated distribution map ofsusceptibility property of the object in the current iteration mayinclude generating the estimated field map in the current iteration byperforming a partial forward dipole kernel function on the updateddistribution map of susceptibility property of the object in the currentiteration.

In some embodiments, the updating the distribution map of susceptibilityproperty of the object in the current iteration based on the errorlimiting information updated in the current iteration may includegenerating the updated distribution map of susceptibility property ofthe object in the current iteration by subtracting an updated artifactestimation map in the current iteration from the distribution map ofsusceptibility property of the object in the preceding iteration; orgenerating the updated distribution map of susceptibility property ofthe object in the current iteration by using a regularization termrelating to an updated artifact estimation map in the iteration process.

In some embodiments, the error limiting information may include an apriori map and an artifact estimation map. The updating the distributionmap of susceptibility property of the object in the current iterationbased on the error limiting information updated in the current iterationmay include generating the updated distribution map of susceptibilityproperty of the object in the current iteration based on the a priorimap in the current iteration, the artifact estimation map in the currentiteration and the preliminary field map.

In some embodiments, the iteration process may be performed based on atleast one of a conjugate gradient, a preconditioned conjugate gradient,quadratic minimization, or split-Bregman.

In some embodiments, the error limiting information may include an apriori map relating to a spatial smoothness, a spatial gradient, edgeinformation, or a weighted property mask of the distribution map ofsusceptibility property of the object.

According to another aspect of the present disclosure, a system isprovided. The system may include at least one storage and at least oneprocessor configured to communicate with the at least one storage. Theat least one storage may include a set of instructions or programs. Whenthe at least one processor executes the set of instructions or programs,the at least one processor may be directed to perform one or more of thefollowing operations. The at least one processor may be directed toobtain a phase diagram corresponding to a magnetic resonance (MR) signalof the object, and determine a preliminary field map based on the phasediagram. The at least one processor may be further directed to obtainpreliminary error limiting information associated with the preliminaryfield map, and determine, based on the preliminary field map and thepreliminary error limiting information, a preliminary distribution mapof susceptibility property of the object. The at least one processor maybe further directed to perform an iteration process including at leastone iteration. In each of the at least one iteration, the at least oneprocessor may be further directed to update error limiting informationin a current iteration based on a distribution map of susceptibilityproperty of the object updated in a preceding iteration; update thedistribution map of susceptibility property of the object in the currentiteration based on the error limiting information updated in the currentiteration; and perform a comparison based on the updated distributionmap in the current iteration and the preliminary field map. The at leastone processor may be further directed to designate the updateddistribution map of susceptibility property of the object as a targetdistribution map of susceptibility property of the object.

According to another aspect of the present disclosure, a computerreadable medium is provided. The computer readable medium may includeexecutable instructions or programs. When executed by at least oneprocessor, the executable instructions or programs may cause the atleast one processor to effectuate a method. The method may include oneor more of the following operations. The method may include obtaining aphase diagram corresponding to a magnetic resonance (MR) signal of theobject and determining a preliminary field map based on the phasediagram. The method may further include obtaining preliminary errorlimiting information associated with the preliminary field map anddetermining, based on the preliminary field map and the preliminaryerror limiting information, a preliminary distribution map ofsusceptibility property of the object. The method may include performingan iteration process including at least one iteration. Each of the atleast one iteration may include updating error limiting information in acurrent iteration based on a distribution map of susceptibility propertyof the object updated in a preceding iteration; updating thedistribution map of susceptibility property of the object in the currentiteration based on the error limiting information updated in the currentiteration; and performing a comparison based on the updated distributionmap in the current iteration and the preliminary field map. The methodmay further include designating the updated distribution map ofsusceptibility property of the object as a target distribution map ofsusceptibility property of the object.

Additional features will be set forth in part in the description whichfollows, and in part will become apparent to those skilled in the artupon examination of the following and the accompanying drawings or maybe learned by production or operation of the examples. The features ofthe present disclosure may be realized and attained by practice or useof various aspects of the methodologies, instrumentalities andcombinations set forth in the detailed examples discussed below.

BRIEF DESCRIPTION OF THE DRAWINGS

The present disclosure is further described in terms of exemplaryembodiments. These exemplary embodiments are described in detail withreference to the drawings. These embodiments are non-limiting exemplaryembodiments, in which like reference numerals represent similarstructures throughout the several views of the drawings, and wherein:

FIG. 1 is a schematic diagram illustrating an exemplary imaging systemaccording to some embodiments of the present disclosure;

FIG. 2 is a schematic diagram illustrating exemplary hardware andsoftware components of a computing device according to some embodimentsof the present disclosure;

FIG. 3 is a schematic diagram illustrating exemplary hardware and/orsoftware components of an exemplary mobile device that is configured toimplement a specific system disclosed in the present disclosure;

FIG. 4 is a schematic diagram illustrating an exemplary scanner and anexemplary processing device according to some embodiments of the presentdisclosure;

FIG. 5 is a schematic diagram illustrating an exemplary processingmodule according to some embodiments of the present disclosure;

FIG. 6 is a flowchart of an exemplary process for determining a targetdistribution map of susceptibility property of an object according tosome embodiments of the present disclosure;

FIG. 7 is a flowchart of an exemplary process for determining apreliminary field map and/or a preliminary a priori map according tosome embodiments of the present disclosure;

FIG. 8 is a flowchart of an exemplary process for iteratively updatingthe distribution map of susceptibility property and/or error limitinginformation according to some embodiments of the present disclosure;

FIG. 9 is a flowchart of an exemplary process for updating the errorlimiting information according to some embodiments of the presentdisclosure;

FIG. 10 is a flowchart of an exemplary process for updating thedistribution map of susceptibility property according to someembodiments of the present disclosure;

FIG. 11 is a flowchart of an exemplary process for determining whether atermination condition relating to the iteration process is satisfiedaccording to some embodiments of the present disclosure;

FIG. 12 is a schematic diagram illustrating an exemplary idealdistribution map of susceptibility property, an exemplary field map, andan exemplary reconstructed distribution map of susceptibility propertywith artifacts according to some embodiments of the present disclosure;

FIG. 13 is a schematic diagram illustrating exemplary schemes ofiteratively updating a priori information according to some embodimentsof the present disclosure;

FIG. 14 is a schematic diagram illustrating exemplary computer simulateddemonstrations of artifact of various amount in the distribution map ofsusceptibility property according to some embodiments of the presentdisclosure; and

FIG. 15 is a schematic diagram illustrating exemplary human brain imagesdemonstrating artifacts of various amount in the distribution maps ofsusceptibility property according to some embodiments of the presentdisclosure.

DETAILED DESCRIPTION

The following description is presented to enable any person skilled inthe art to make and use the present disclosure, and is provided in thecontext of a particular application and its requirements. Variousmodifications to the disclosed embodiments will be readily apparent tothose skilled in the art, and the general principles defined herein maybe applied to other embodiments and applications without departing fromthe spirit and scope of the present disclosure. Thus, the presentdisclosure is not limited to the embodiments shown, but is to beaccorded the widest scope consistent with the claims.

It will be understood that the term “system,” “engine,” “unit,”“module,” and/or “block” used herein are one method to distinguishdifferent components, elements, parts, section or assembly of differentlevel in ascending order. However, the terms may be displaced by anotherexpression if they may achieve the same purpose.

Generally, the word “module,” “unit,” or “block,” as used herein, refersto logic embodied in hardware or firmware, or to a collection ofsoftware instructions. A module, a unit, or a block described herein maybe implemented as software and/or hardware and may be stored in any typeof non-transitory computer-readable medium or another storage device. Insome embodiments, a software module/unit/block may be compiled andlinked into an executable program. It will be appreciated that softwaremodules can be callable from other modules/units/blocks or fromthemselves, and/or may be invoked in response to detected events orinterrupts. Software modules/units/blocks configured for execution oncomputing devices (e.g., CPU 220 as illustrated in FIG. 2) may beprovided on a computer-readable medium, such as a compact disc, adigital video disc, a flash drive, a magnetic disc, or any othertangible medium, or as a digital download (and can be originally storedin a compressed or installable format that needs installation,decompression, or decryption prior to execution). Such software code maybe stored, partially or fully, on a storage device of the executingcomputing device, for execution by the computing device. Softwareinstructions may be embedded in firmware, such as an ErasableProgrammable Read Only Memory (EPROM). It will be further appreciatedthat hardware modules/units/blocks may be included in connected logiccomponents, such as gates and flip-flops, and/or can be included ofprogrammable units, such as programmable gate arrays or processors. Themodules/units/blocks or computing device functionality described hereinmay be implemented as software modules/units/blocks, but may berepresented in hardware or firmware. In general, themodules/units/blocks described herein refer to logicalmodules/units/blocks that may be combined with othermodules/units/blocks or divided into sub-modules/sub-units/sub-blocksdespite their physical organization or storage. The description may beapplicable to a system, an engine, or a portion thereof.

The terminology used herein is for the purpose of describing particularexample embodiments only and is not intended to be limiting. As usedherein, the singular forms “a,” “an,” and “the” may be intended toinclude the plural forms as well, unless the context clearly indicatesotherwise. It will be further understood that the terms “comprise,”“comprises,” and/or “comprising,” “include,” “includes,” and/or“including,” when used in this specification, specify the presence ofstated features, integers, steps, operations, elements, and/orcomponents, but do not preclude the presence or addition of one or moreother features, integers, steps, operations, elements, components,and/or groups thereof. It should be noted that the terms “magneticsusceptibility map,” “susceptibility map,” “susceptibilitydistribution,” “susceptibility distribution map,” and “distribution mapof susceptibility property” may be used interchangeably in the presentdisclosure to refer to a 3D image or image data of the susceptibilityproperty of a scanned object. It should also be noted that the terms “apriori field map” and “estimated field map” may be used interchangeablyin the present disclosure to refer to a 3D image or image data of themagnetic field around an object.

These and other features, and characteristics of the present disclosure,as well as the methods of operation and functions of the relatedelements of structure and the combination of parts and economies ofmanufacture, may become more apparent upon consideration of thefollowing description with reference to the accompanying drawings, allof which form a part of this disclosure. It is to be expresslyunderstood, however, that the drawings are for the purpose ofillustration and description only and are not intended to limit thescope of the present disclosure. It is understood that the drawings arenot to scale.

The flowcharts used in the present disclosure illustrate operations thatsystems implement according to some embodiments in the presentdisclosure. It is to be expressly understood, the operations of theflowchart may be implemented not in order. Conversely, the operationsmay be implemented in inverted order, or simultaneously. Moreover, oneor more other operations may be added to the flowcharts. One or moreoperations may be removed from the flowcharts.

The present disclosure relates to a method and system for determining adistribution map of susceptibility property of an object. For example, amagnetic resonance (MR) signal of the object may be obtained. The MRsignal may include intensity information and frequency information. Amagnitude diagram may be generated based on the intensity information,and a phase diagram may be generated based on the frequency information.A preliminary field map may be generated based on the phase diagram. Apreliminary error limiting information may be obtained. The preliminaryerror limiting information may include preliminary a priori informationand/or preliminary artifact estimation information. A preliminarydistribution map of susceptibility property of the object may bedetermined based on the preliminary field map and the preliminary errorlimiting information. An iterative process may be performed. Theiterative process may include a plurality of iterations. In eachiteration, the error limiting information and the distribution map ofsusceptibility property may be updated. After the updating of the errorlimiting information and the distribution map of susceptibility propertyin each iteration, a difference field map may be constructed based onthe updated distribution map of susceptibility property and thepreliminary field map. If an evaluation value related to the differencefield map is less than a threshold, the iterative process may terminate,and the updated distribution map of susceptibility property may bedesignated as a target distribution map of susceptibility property ofthe object.

If the error limiting information includes a prior information, in theiterative process for updating the distribution map of susceptibilityproperty, the a prior information can be dynamically extracted andupdated based on the intermediate distribution map of susceptibilityproperty, and the updated a prior information may be applied to theupdating of the distribution map of susceptibility property, therebymaking the a prior information approach the real spatial characteristicsof the susceptibility distribution, and accelerating the convergencespeed of the iterative process and improve the accuracy of the results(i.e., the distribution map of susceptibility property).

If the error limiting information includes artifact estimationinformation, in the iterative process, the artifact estimationinformation is dynamically updated based on intermediate distributionmap of susceptibility property, and the updated artifact estimationinformation may be used as a regularization term and applied to theiterative process (or the updated artifact estimation information may beremoved directly from the distribution map of susceptibility property),thereby reducing the artifact intensity of the distribution map ofsusceptibility property to a negligible level without affecting othercalculation procedures, and ensuring the accuracy of the results (i.e.,the distribution map of susceptibility property).

FIG. 1 is a schematic diagram illustrating an exemplary imaging systemaccording to some embodiments of the present disclosure. As shown inFIG. 1, the imaging system 100 may include a scanner 110, a network 120,one or more terminals 130, a processing device 140, and a storage device150. All the components in the imaging system 100 may be interconnectedvia the network 120.

The scanner 110 may scan an object 114 and generate scanned datarelating to the object 114. The object 114 may be biological ornon-biological. Merely by way of example, the object 114 may include apatient, an organ, a tissue, a specimen, a man-made object, a phantom,etc. In some embodiments, the scanner 110 may be a single-modalitymedical imaging device, such as, an MRI device, a PET device, a SPECTdevice, a CT device, or the like, or a multi-modality medical imagingdevice, such as, a PET-MRI device, a SPECT-MRI device, or a PET-CTdevice. Taking an MRI scanner (e.g., the MRI device, the PER-MRI device,the SPECT-MRI device) as an example, the scanner 110 may include agantry 111 and a scanning bed 113. The scanner 110 may also include amain magnet, a gradient magnet, a pulse controller, a gradient signalgenerator, a radiofrequency (RF) pulse generator, an RF signal receiver(or collectively referred to as an RF signal transceiver), and coils(not shown in FIG. 1). The scanner 110 may obtain information related tothe object 114 by scanning the object 114. More particularly, certainatomic nuclei (e.g., hydrogen-1, carbon-13, oxygen-17, etc.) may absorband emit RF energy when being placed in the main magnetic and thegradient magnetic fields. The emitted RF energy exists as RF signals andis received by the scanner 110. The scanner 110 may obtain informationfor reconstructing a distribution of the atomic nuclei inside the object14 based on the RF signals as an MRI image.

The main magnet may be placed in the gantry 111 of the scanner 110. Themain magnet may be configured to generate a (substantially) uniform mainmagnetic field. The strength of the main magnetic field may be, e.g.,0.2 Tesla, 0.5 Tesla, 1.0 Tesla, 1.5 Tesla, 3.0 Tesla, etc. In someembodiments, the main magnet may be a superconducting coil.Alternatively, the main magnet may be a permanent magnet.

The scanning bed 113 may support the object 114 during a scan. Forexample, the object 114 may be supported and/or delivered to thedetecting region 112 of the gantry 111 by the scanning bed 113. Thedetecting region 112 may be a region that the magnetic fielddistribution of the main magnetic field is relatively uniform, and towhich, the RF signal is transmitted.

Merely by way of example, a spatial coordinate system (i.e., acoordinate system of the scanner 110) may be described with reference tothe gantry 111 of the scanner 110. For example, a Z axis may be thedirection along the axis of the gantry, an X axis and a Y axis may bethe directions perpendicular to the Z axis. As illustrated in FIG. 1,the X axis may be along the horizontal direction, and the Y axis may bealong the vertical direction. The long-axis direction of the object 114may coincide with the direction of the Z axis, and the scanning bed 113may move in the direction of the Z axis.

The pulse controller may be configured to control the generation of RFsignals. For example, the pulse controller may control the RF pulsegenerator and the gradient signal generator. In some embodiments, thepulse controller may control the RF pulse generator to generate an RFpulse. The RF pulse may be amplified by an amplifier. In someembodiments, the pulse controller may control the gradient signalgenerator to generate gradient signals. The pulse controller may receiveinformation from or send information to the scanner 110, the processingdevice 140, and/or a terminal 130. According to some embodiments of thepresent disclosure, the pulse controller may receive a command from auser via the terminal 130 or a display (e.g., display 450) of theprocessing device 140 and control components of the scanner 110 (e.g.,the RF pulse generator) accordingly to start a scan.

The RF pulse generator may generate an RF pulse. In some embodiments,the RF pulse generator may generate the RF pulse based on an instructionfrom the pulse controller. The RF pulse may be amplified by anamplifier. A switch may be configured to control the emission of theamplified RF pulse. For example, the amplified RF pulse may be emittedby the coils when the switch is on. The emitted RF pulse may excite theatomic nuclei in the object 114. The object 114 may generate acorresponding RF signal when the RF excitation is removed. The coils maytransmit RF signals to or receive RF signals from the object 114. Forexample, the coils may transmit the amplified RF pulse to the object114. As another example, the coils may receive the RF signals emittedfrom the object 114. In some embodiments, the coils may include aplurality of RF receiving channels. The plurality of RF receivingchannels may transmit RF signals received from the object 114 to the RFsignal receiver. The RF signal receiver may be configured to receive RFsignals. The RF signal receiver may receive RF signals from the coils.In some embodiments, the coils may include a volume coil configured totransmit RF signals to or receive RF signals from the entire body of theobject 114 and/or a local coil configured to transmit RF signals to orreceive RF signals from a portion of the object 114.

The gradient signal generator may generate gradient signals. In someembodiments, the gradient signal generator may generate gradient signalsbased on an instruction from the pulse controller. The gradient signalsmay include three orthogonal signals. In some embodiments, the threeorthogonal signals may be a first gradient signal along the X direction,a second gradient signal along the Y direction, and a third gradientsignal along the Z direction. The gradient signals may facilitate tolocate the atomic nuclei.

The gradient magnet may be configured to spatially encode RF signals(e.g., an RF pulse generated by the RF pulse generator). The gradientmagnet may generate a magnetic field with a strength less than that ofthe main magnetic field. For example, the gradient magnet may generate agradient magnetic field in the detecting region 112.

The network 120 may include any suitable network that can facilitate theexchange of information and/or data for the imaging system 100. In someembodiments, one or more components of the imaging system 100 (e.g., thescanner 110, the terminal 130, the processing device 140, the storagedevice 150, etc.) may communicate information and/or data with one ormore other components of the imaging system 100 via the network 120. Forexample, the processing device 140 may obtain image data from thescanner 110 via the network 120. As another example, the processingdevice 140 may obtain user instructions from the terminal 130 via thenetwork 120. The network 120 may be and/or include a public network(e.g., the Internet), a private network (e.g., a local area network(LAN), a wide area network (WAN)), etc.), a wired network (e.g., anEthernet network), a wireless network (e.g., an 802.11 network, a Wi-Finetwork, etc.), a cellular network (e.g., a Long Term Evolution (LTE)network), a frame relay network, a virtual private network (“VPN”), asatellite network, a telephone network, routers, hubs, switches, servercomputers, and/or any combination thereof. Merely by way of example, thenetwork 120 may include a cable network, a wireline network, afiber-optic network, a telecommunications network, an intranet, awireless local area network (WLAN), a metropolitan area network (MAN), apublic telephone switched network (PSTN), a Bluetooth™ network, aZigBee™ network, a near field communication (NFC) network, or the like,or any combination thereof. In some embodiments, the network 120 mayinclude one or more network access points. For example, the network 120may include wired and/or wireless network access points such as basestations and/or internet exchange points through which one or morecomponents of the imaging system 100 may be connected to the network 120to exchange data and/or information.

The terminal(s) 130 may include a mobile device 131, a tablet computer132, a laptop computer 133, or the like, or any combination thereof. Insome embodiments, the mobile device 130-1 may include a smart homedevice, a wearable device, a mobile device, a virtual reality device, anaugmented reality device, or the like, or any combination thereof. Insome embodiments, the smart home device may include a smart lightingdevice, a control device of an intelligent electrical apparatus, a smartmonitoring device, a smart television, a smart video camera, aninterphone, or the like, or any combination thereof. In someembodiments, the wearable device may include a bracelet, a footgear,eyeglasses, a helmet, a watch, clothing, a backpack, a smart accessory,or the like, or any combination thereof. In some embodiments, the mobiledevice may include a mobile phone, a personal digital assistant (PDA), agaming device, a navigation device, a point of sale (POS) device, alaptop, a tablet computer, a desktop, or the like, or any combinationthereof. In some embodiments, the virtual reality device and/or theaugmented reality device may include a virtual reality helmet, virtualreality glasses, a virtual reality patch, an augmented reality helmet,augmented reality glasses, an augmented reality patch, or the like, orany combination thereof. For example, the virtual reality device and/orthe augmented reality device may include a Google Glass™, an OculusRift™, a Hololens™, a Gear VR™, etc. In some embodiments, theterminal(s) 130 may be part of the processing device 140.

The processing device 140 may process data and/or information obtainedfrom the scanner 110, the terminal 130, and/or the storage device 150.In some embodiments, the processing device 140 may be a single server ora server group. The server group may be centralized or distributed. Insome embodiments, the processing device 140 may be local or remote. Forexample, the processing device 140 may access information and/or datastored in the scanner 110, the terminal 130, and/or the storage device150 via the network 120. As another example, the processing device 140may be directly connected to the scanner 110, the terminal 130 and/orthe storage device 150 to access stored information and/or data. In someembodiments, the processing device 140 may be implemented on a cloudplatform. Merely by way of example, the cloud platform may include aprivate cloud, a public cloud, a hybrid cloud, a community cloud, adistributed cloud, an inter-cloud, a multi-cloud, or the like, or anycombination thereof. In some embodiments, the processing device 140 maybe implemented by a computing device 200 having one or more componentsas illustrated in FIG. 2. The processing device 140 may be configured toprocess the MR signals received from the RF signal receiver and/orgenerate one or more MR images based on the MR signals. The processingdevice 140 may process the received MR signals to generate and/or updatea magnitude diagram, a phase diagram, a susceptibility map, an R2* (T2*)map, a susceptibility weighted imaging (SWI) map, or the like, or acombination thereof. In some embodiments, the processing device 140 mayperform an iterative process for updating the distribution map ofsusceptibility property and error limiting information.

The storage device 150 may store data, instructions, and/or any otherinformation. In some embodiments, the storage device 150 may store dataobtained from the terminal 130 and/or the processing device 140. In someembodiments, the storage device 150 may store data and/or instructionsthat the processing device 140 may execute or use to perform exemplarymethods described in the present disclosure. In some embodiments, thestorage device 150 may include mass storage, removable storage, avolatile read-and-write memory, a read-only memory (ROM), or the like,or any combination thereof. Exemplary mass storage may include amagnetic disk, an optical disk, a solid-state drive, etc. Exemplaryremovable storage may include a flash drive, a floppy disk, an opticaldisk, a memory card, a zip disk, a magnetic tape, etc. Exemplaryvolatile read-and-write memories may include a random access memory(RAM). Exemplary RAM may include a dynamic RAM (DRAM), a double daterate synchronous dynamic RAM (DDR SDRAM), a static RAM (SRAM), athyristor RAM (T-RAM), and a zero-capacitor RAM (Z-RAM), etc. ExemplaryROM may include a mask ROM (MROM), a programmable ROM (PROM), anerasable programmable ROM (EPROM), an electrically erasable programmableROM (EEPROM), a compact disk ROM (CD-ROM), and a digital versatile diskROM, etc. In some embodiments, the storage device 150 may be implementedon a cloud platform. Merely by way of example, the cloud platform mayinclude a private cloud, a public cloud, a hybrid cloud, a communitycloud, a distributed cloud, an inter-cloud, a multi-cloud, or the like,or any combination thereof.

In some embodiments, the storage device 150 may be connected to thenetwork 120 to communicate with one or more other components in theimaging system 100 (e.g., the processing device 140, the terminal 130,etc.). One or more components of the imaging system 100 may access thedata or instructions stored in the storage device 150 via the network120. In some embodiments, the storage device 150 may be directlyconnected to or communicate with one or more other components of theimaging system 100 (e.g., the processing device 140, the terminal 130,etc.). In some embodiments, the storage device 150 may be part of theprocessing device 140.

FIG. 2 is a schematic diagram illustrating exemplary hardware andsoftware components of a computing device according to some embodimentsof the present disclosure. The computing device 200 may be a generalpurpose computer or a special purpose computer; both may be used toimplement an imaging system 100 of the present disclosure. In someembodiments, the processing device 140 may be implemented on thecomputing device 200, via its hardware, software program, firmware, or acombination thereof. Although only one such computer is shown, forconvenience, the computer functions as described herein may beimplemented in a distributed manner on a number of similar platforms, todistribute the processing load.

The computing device 200, for example, may include COM ports 250connected to and from a network connected thereto to facilitate datacommunications. The computing device 200 may also include a centralprocessing unit (CPU) 220, in the form of one or more processors, forexecuting program instructions. The exemplary computer platform mayinclude an internal communication bus 210, program storage and datastorage of different forms, for example, a disk 270, and a read-onlymemory (ROM) 230, or a random access memory (RAM) 240, for various datafiles to be processed and/or transmitted by the computer. The exemplarycomputer platform may also include program instructions stored in theROM 230, RAM 240, and/or another type of non-transitory storage mediumto be executed by the CPU 220. The methods and/or processes of thepresent disclosure may be implemented as the program instructions. Thecomputing device 200 also includes an I/O component 260, supportinginput/output between the computer and other components therein such asuser interface elements 280. The computing device 200 may also receiveprogramming and data via network communications.

The computing device 200 may also include a hard disk controllercommunicated with a hard disk, a keypad/keyboard controller communicatedwith a keypad/keyboard, a serial interface controller communicated witha serial peripheral equipment, a parallel interface controllercommunicated with a parallel peripheral equipment, a display controllercommunicated with a display, or the like, or any combination thereof.

Merely for illustration, only one CPU and/or processor is described inthe computing device 200. However, it should be noted that the computingdevice 200 in the present disclosure may also include multiple CPUsand/or processors, thus operations and/or method steps that areperformed by one CPU and/or processor as described in the presentdisclosure may also be jointly or separately performed by the multipleCPUs and/or processors. For example, if in the present disclosure theCPU and/or processor of the computing device 200 executes both operationA and operation B (e.g., the updating of distribution map ofsusceptibility property of the object and the updating of error limitinginformation), it should be understood that operation A and operation Bmay also be performed by two or more different CPUs and/or processorsjointly or separately in the computing device 200 (e.g., a firstprocessor updates the distribution map of susceptibility property of theobject and a second processor updates the error limiting information, orthe first and second processors jointly update the distribution map ofsusceptibility property of the object and the error limitinginformation).

FIG. 3 is a schematic diagram illustrating exemplary hardware and/orsoftware components of an exemplary mobile device that is configured toimplement a specific system disclosed in the present disclosure. Asillustrated in FIG. 3, the mobile device 300 may include a communicationunit 310, a display 320, a graphics processing unit (GPU) 330, a CPU340, an I/O 350, a storage 390, and a memory 360. In some embodiments,any other suitable component, including but not limited to a system busor a controller (not shown), may also be included in the mobile device300. In some embodiments, a mobile operating system 370 (e.g., IOS™,Android™, Windows Phone™, etc.) and one or more applications 380 may beloaded into the memory 360 from the storage 390 in order to be executedby the CPU 340. The applications 380 may include a browser or any othersuitable mobile apps for receiving and rendering information relating toimage processing or other information from the processing device 140.User interactions with the information stream may be achieved via theI/O 350 and provided to the processing device 140 and/or othercomponents of the imaging system 100 via the network 120. In someembodiments, a user may input parameters to the imaging system 100, viathe mobile device 300.

In order to implement various modules, units and their functionsdescribed above, a computer hardware platform may be used as hardwareplatforms of one or more elements (e.g., the processing device 140and/or other components of the imaging system 100 described in FIG. 1).Since these hardware elements, operating systems and program languagesare common; it may be assumed that persons skilled in the art may befamiliar with these techniques and they may be able to provideinformation required in the imaging according to the techniquesdescribed in the present disclosure. A computer with the user interfacemay be used as a personal computer (PC), or other types of workstationsor terminal devices. After being properly programmed, a computer withthe user interface may be used as a server. It may be considered thatthose skilled in the art may also be familiar with such structures,programs, or general operations of this type of computer device. Thus,extra explanations are not described for the Figures.

FIG. 4 is a schematic diagram illustrating an exemplary scanner and anexemplary processing device according to some embodiments of the presentdisclosure. As shown in FIG. 4, the scanner 110 may include a magneticfield generating module 410 and a radiofrequency signal transceivingmodule 420. The processing device 140 may include a controlling module430, a processing module 440, and a display 450. In some embodiments,the scanner 110 may communicate with the processing device 140 viacables. Alternatively or additionally, the scanner 110 may communicatewith the processing device 140 via the network 120. In some embodiments,the scanner 110 and the processing device 140 may be integrated into asingle component that performs either or both of their functions.

In some embodiments, the scanner 110 may perform a scan on an object(e.g., the object 114) to obtain/acquire magnetic resonance (MR) signalsof the object's target site (or region of interest (ROI)). In someembodiments, the scan may be a normal MRI scan on the object (e.g., withnormal settings and/or parameters) for generating an MR image of theobject or a pre-scan of a phantom for calibrating the imaging system100.

The magnetic field generating module 410 may include a main magneticfield generator and/or a gradient magnetic field generator. The mainmagnetic field generator may generate a static magnetic field B₀ duringthe scan. The main magnetic field generator may be of various typesincluding, for example, a permanent magnet, a superconductive magnet, aresistive electromagnet, or the like, or a combination thereof. Thegradient magnetic field generator may generate magnetic field gradientsGx, Gy, and Gz in the directions of “X,” “Y,” and “Z,” respectively. Asused herein, the X, Y, and Z directions may represent the directionsalong the X, Y, and Z axes in a coordinate system associated with thescanner 110. Similar to FIG. 1, the X and Z axes may be in a horizontalplane, the X and Y axes may be in a vertical plane, and the Z axis maybe along the axis of rotation of the gantry. In some embodiments, the Xaxis, the Y axis, and the Z axis may be specified by the gradientmagnetic field generator (e.g., a gradient coil in a gradient magneticfield generator). The gradient magnetic field generator may encodeand/or read spatial information (e.g., X, Y, Z coordinates) of an objector a portion thereof located within the scanner 110. In someembodiments, the gradient magnetic field generator may generate one ormore gradient magnetic fields, each of which is in a particulardirection during the scan. Merely by way of example, the gradientmagnetic field generator may generate a first gradient magnetic fieldalong the X axis, a second gradient magnetic field along the Y axis, anda third gradient magnetic field along the Z axis. In some embodiments,the gradient magnetic fields along the X axis, the Y axis, and/or the Zaxis may correspond to different encoding/reading directions in k-space(e.g., the direction of the Kx axis, the direction of the Ky axis, thedirection of the Kz axis, or any other directions).

The radiofrequency signal transceiving module 420 may include an RFtransmitting coil and/or an RF receiving coil. The RF transmitting coilmay transmit RF signals to the object and the RF receiving coil mayreceive RF signals from the object. In some embodiments, thefunctionality, size, type, geometry, location, amount, and/or magnitudeof the magnetic field generating module 410 and/or the radiofrequencysignal transceiving module 420 may be changed according to one or morespecific requirements in different scenarios. For example, RF coils maybe classified into volume coils and local coils depending on differencesin, e.g., their function, size, or the like, or a combination thereof.Exemplary volume coils may include a birdcage coil, a transverseelectromagnetic coil, a surface coil, a saddle coil, or the like, or acombination thereof. Exemplary local coils may include a birdcage coil,a solenoid coil, a saddle coil, a flexible coil, or the like, or acombination thereof. In some embodiments, the magnetic field generatingmodule 410 and the radiofrequency signal transceiving module 420 may beconfigured to surround the object to form a tunnel-type MRI scanner(e.g., a closed-cell MRI scanner, such as the scanner 110 shown inFIG. 1) or an open MRI scanner.

The controlling module 430 may control the magnetic field generatingmodule 410 and/or the radiofrequency signal transceiving module 420 ofthe scanner 110, the processing module 440, and/or the display 450 ofthe processing device 140, and/or other components of the imaging system100. For example, the controlling module 430 may control the strengthand direction of the gradient magnetic fields generated by the gradientmagnetic field generator. As another example, the controlling module 430may receive or transmit information from/to the scanner 110, theprocessing module 440, and/or the display 450. In some embodiments, thecontrolling module 430 may receive a command (e.g., provided by a user)from the display 450, and operate the magnetic field generating module410 and/or the radiofrequency signal transceiving module 420 accordingto the received command to obtain an image of the object. Merely by wayof example, the command may relate to an imaging application (e.g.,selected by the user), which may include but is not limited tosusceptibility weighted imaging (SWI), or quantitative susceptibilitymapping (QSM), or the like.

The processing module 440 may be configured to process various kinds ofinformation received from different modules. In some embodiments, theprocessing module 440 may include a central processing unit (CPU), anapplication specific integrated circuit (ASIC), an application-specificinstruction-set processor (ASIP), a graphics processing unit (GPU), aphysics processing unit (PPU), a microcontroller, a digital signalprocessor (DSP), a field programmable gate array (FPGA), an advancedRISC machine (ARM), or the like, or a combination thereof. Theprocessing module 440 may process the MR signals received from theradiofrequency signal transceiving module 420 and generate one or moreMR images based on the MR signals. The processing module 440 may deliverthe images to the display 450. In some embodiments, the processingmodule 440 may process data inputted by the user or operator via thedisplay 450 and transform the data into specific commands and transmitthe commands to the controlling module 430. In some embodiments, theprocessing module 440 may perform the processes as illustrated in FIG. 6through FIG. 11 to process the received magnetic resonance signals togenerate and/or update a magnitude diagram, a phase diagram, asusceptibility map, an R2* (T2*) map, an SWI map, or the like, or acombination thereof.

The display 450 may be configured to display output information. In someembodiments, an input device may be configured to allow a user toprovide input. In some embodiments, the input device may be implementedon or operably connected to the display 450. For instance, the display450 may be a touch screen that allows input and display. The inputinformation and/or output information may include programs, software,algorithms, data, texts, numbers, images, sounds, or the like, or anycombination thereof. For example, the user or operator may enter initialMR parameters or settings to initiate a scan. In some embodiments, thecontrolling module 430, the processing module 440, and/or the display450 may be integrated into an image generator (e.g., an independentimage generator, the terminals 130). The user may set parameters for theMR scan, control the imaging protocol, view the images produced by theimage generator on a display thereof. In some embodiments, the display450 may be a CRT display or an LCD. The display 450 may display an MRsignal, a magnitude diagram, a phase diagram, a susceptibility map, anR2* (T2*) map, an SWI map, or the like, or a combination thereof.

FIG. 5 is a schematic diagram illustrating an exemplary processingmodule according to some embodiments of the present disclosure. As shownin FIG. 5, the processing module 440 may include an acquisition unit510, an updating unit 520, and a processing unit 530.

The acquisition unit 510 may be configured to obtain a preliminarydistribution map of susceptibility property of the object based on thepreliminary field map and the preliminary error limiting information.Specifically, the preliminary distribution map of susceptibilityproperty of the object may be obtained based on the preliminary fieldmap and at least one of the preliminary a priori information or thepreliminary artifact estimation information. In some embodiments, theacquisition unit 510 may receive the MR signals acquired by the scanner110 and determine the preliminary field map based on the MR signals ofthe target site of the object. In some embodiments, the acquisition unit510 may obtain at least one of a magnitude diagram, a phase diagram, anSWI map, or an R2* (T2*) map of the object. The acquisition unit 510 maydetermine the preliminary a priori information based on the at least oneof the magnitude diagram, the phase diagram, the SWI map, or the R2*(T2*) map. The preliminary a priori information may include a spatialsmoothness of the susceptibility distribution, a spatial gradient of thesusceptibility distribution, edge information of the susceptibilitydistribution, a weighted property mask, or the like, or any combinationthereof.

The updating unit 520 may be configured to perform an iterative processfor updating the distribution map of susceptibility property and/or theerror limiting information. The iterative process may include aplurality of iterations. In some embodiments, in each iteration, theerror limiting information in the current iteration may be updated basedon the distribution map of susceptibility property in the precedingiteration, and the distribution map of susceptibility property in thecurrent iteration may be updated based on the error limiting informationupdated in the current iteration. The error limiting information mayinclude the a priori information and/or the remnant streaking artifactestimation information corresponding to the preliminary field map. Forexample, the updating unit 520 may, in each iteration, update the apriori information, and determine the distribution map of susceptibilityproperty based on the updated a priori information. As another example,the updating unit 520 may, in each iteration, update the remnantstreaking artifact estimation information corresponding to the fieldmap, and determine the distribution map of susceptibility property basedon the updated remnant radiation streaking artifacts corresponding tothe field map. As a further example, the updating unit 520 may, in eachiteration, update both the a priori information and the remnantstreaking artifact estimation information, either in sequence orsimultaneously, and determine the distribution map of susceptibilityproperty based on the updated a priori information and the updatedremnant streaking artifact estimation information.

The processing unit 530 may obtain the updated distribution map ofsusceptibility property from the updating unit 520. The iterativeprocess may terminate when a termination condition is satisfied. Forinstance, the termination condition may be that the evaluation valuebetween the preliminary field map and a field map corresponding to theupdated distribution map of susceptibility property (e.g., a field mapgenerated by performing a forward dipole kernel function on the updateddistribution map of susceptibility property) is less than a thresholdwhen the iterative process terminates. The processing unit 530 mayoutput the updated distribution map of susceptibility property as atarget distribution map of susceptibility property of the object.

FIG. 6 is a flowchart of an exemplary process for determining a targetdistribution map of susceptibility property of an object according tosome embodiments of the present disclosure. The process 600 may beexecuted by the processing device 140. For example, the process 600 maybe implemented as a set of instructions (e.g., an application) stored inthe storage, e.g., ROM 230, RAM 240, the storage device 150, the storage390, a storage device external to and accessible by the imaging system100. The processing device 140, the CPU 220, the CPU 340, the processingmodule 440, and/or the processing unit 530 may execute the set ofinstructions, and when executing the instructions, it may be configuredto perform the process 600.

In 602, the acquisition unit 510 may receive a magnetic resonance (MR)signal from an object (e.g., the object 114). In some embodiments, theobject may include a patient, a tissue or organ of a patient, asubstance, a specimen, or the like, or any combination thereof. Forexample, the object may include the head, the chest, the lung, thepleura, the mediastinum, the abdomen, the large intestine, the smallintestine, the bladder, the gallbladder, the chest, the pelvic cavity,the backbone, the limb(s), the skeleton, the blood vessel(s), etc., of apatient. In some embodiments, the object or a region of interest (ROI,also referred to as target site) of the object may be excited by an RFsignal. When the object or an ROI of the object relaxes, a spin signalor spin echo may be emitted by the object or the ROI of the object. Thespin signal or spin echo may be received by the imaging system 100. Thespin signal may be phase encoded to form a magnetic resonance (MR)signal.

In 604, the processing unit 530 may determine a preliminary field map ofthe object based on the MR signal. In some embodiments, phaseinformation of the MR signal may be extracted to generate thepreliminary field map of the object. More descriptions regarding thegeneration of the preliminary field map of the object may be found in,e.g., operations 704 and 706 in FIG. 7 and the descriptions thereof.

In 606, the acquisition unit 510 may obtain preliminary error limitinginformation of the object. The preliminary error limiting informationmay include but not limited to a priori information and artifactestimation information.

The a priori information may provide features on the spatialdistribution of susceptibility of the object, such as the continuity ofdistribution within the same tissue, the boundaries between differenttissues, etc. For example, the a priori information may include spatialsmoothness of the susceptibility distribution of the object, a spatialgradient of the susceptibility distribution of the object, edgeinformation of the susceptibility distribution of the object, a weightedproperty mask, or the like, or any combination thereof, corresponding tothe target site of the object. The a priori information may be obtainedfrom a magnitude map, a phase map, a magnetic susceptibility map, an R2*(T2*) map, or a susceptibility weighted imaging (SWI) map of the object.More descriptions regarding the generation of the preliminary a prioriinformation may be found in, e.g., operations 708, 710 and 712 in FIG. 7and the descriptions thereof.

The artifact estimation information may include radial streakingartifacts. Due to the ill-posed nature of the inverse dipole problem(e.g., singularity values in the inverse dipole kernel functions),computationally relevant radial streaking artifacts are generated duringQSM. Removing the remnant streaking artifacts in the susceptibility mapmay effectively reduce the intensity of the artifacts and improve theaccuracy of the susceptibility map. In some embodiments, preliminaryartifact estimation information may be generated based on an estimation.In a case that a preliminary distribution map of susceptibility propertyof the object is known (e.g., obtained based on the preliminary fieldmap and the a priori information), the preliminary artifact estimationinformation may be obtained based on a difference map between thepreliminary field map and a field map corresponding to the preliminarydistribution map of susceptibility property (e.g., an estimated fieldmap, a field map generated by performing a forward dipole kernelfunction on the preliminary distribution map of susceptibilityproperty).

In 608, the processing unit 530 may determine a preliminary distributionmap of susceptibility property of the object based on the preliminaryfield map and the preliminary error limiting information. Specifically,the preliminary distribution map of susceptibility property of theobject may be obtained based on the preliminary field map and at leastone of the preliminary a priori information and the preliminary artifactestimation information.

In some embodiments, the preliminary distribution map of susceptibilityproperty of the object may be determined based on Formula (1) asillustrated below:argmin_(χ) ∥F ⁻¹ DFχ−ϕ∥ ₂ ² +α∥P∇χ∥ ₁ +β∥F ⁻¹ D′FΔϕ∥ ₂ ²,  (1)where the first part, ∥F⁻¹DFχ−ϕ∥₂ ², is a basic quadratic minimizationof a conversion function between a field map and a susceptibilitydistribution. In the first part, F and F⁻¹ denote Fourier and inverseFourier operators, respectively, D denotes the forward dipole kernel, χdenotes the susceptibility distribution, and ϕ denotes the fielddistribution (i.e., the field map). In operation 608, χ denotes thepreliminary distribution map of susceptibility property of the object,and ϕ denotes the preliminary field map.

The second part, α∥P∇Δχ∥₁, is the L1 norm of gradients in all threedimensions of the a priori information. In the second part, α denotes aregularization factor that balances data consistency and spatialsmoothness, P denotes a (weighted) property mask that is dynamicallyextracted from the initial or intermediate χ during an iterativeoptimization process, and ∇ denotes a gradient operator. The second partindicates a dynamic updating of the a priori information and/or the useof dynamically updated a priori information. In operation 608, P denotesthe preliminary a priori information (e.g., a preliminary propertymask). In some embodiments, the property mask may be an image includinga plurality of pixels whose values are, e.g., zeros and ones. Forinstance, values of pixels in a smooth region (e.g., a region inside oroutside a type of tissue of the object that have a relatively lowgradient variance) may be zero, while values of pixels near boundaries(e.g., boundaries of different tissues, boundaries of different portionsof a same tissue) having a relatively high gradient variance may be one.For example, a smooth region inside or outside the liver of the objectmay have a homogeneous structure, and the gradient variance of thesmooth region may be relatively low. As another example, the boundary ofthe liver of the object may be relatively sharp to differentiate theregion inside the liver from that outside the liver, the gradients ofthe pixels near the boundary may have a large saltation, and then thegradient variance may be relatively high. In some embodiments, thevalues of pixels may be any value between zero and one. For example, theabove-mentioned pixels inside or outside a type of tissue with arelatively low gradient variance may have any value between zero and one(e.g., 0.01, 0.1, 0.3), and the above-mentioned pixels near a boundarybetween the region inside the tissue and that outside the tissue may beany number between zero and one (e.g., 0.99, 0.9, 0.7). In someembodiments, the values of pixels in the property mask P may be anypositive number (even greater than 1).

The third part, β∥F⁻¹D′FΔϕ∥₂ ², is the L2 norm of the estimation of theremnant streaking artifact. In the third part, β denotes a scalingfactor that adjusts the weighting of the remnant streaking artifactregularization, and D′ denotes a partial forward dipole kernel. D′ mayinclude only the parts in the dipole kernel that include or are close tothe singularities. In some embodiments, the partial forward dipolekernel D′ may be obtained by a threshold calculation based on D (e.g.,D′=D·(|D|<thresh)). Δϕ in the third part denotes the difference (alsoreferred to as a difference field map) between an intermediate (orestimated) field map (calculated from the preliminary distribution mapof susceptibility property or an intermediate χ map using the forwarddipole kernel) and the preliminary field map. The third part indicates adynamic estimation of streaking artifact and mitigation of the streakingartifact, and the use of dynamically updated streaking artifact. Inoperation 608, Δϕ denotes the difference between the preliminary fieldmap and the estimated field map corresponding to the preliminarydistribution map of susceptibility property.

Formula (1) describes an optimization process, e.g., seeking for a valueχ to make the value of the whole Formula (1) minimum or less than athreshold. Formula (1) may be solved by any kinds of numericoptimization techniques, such as a conjugate gradient, a preconditionedconjugate gradient, quadratic minimization, split-Bregman, or the like,or any combination thereof.

It should be noted that Formula (1) is only an exemplary way ofobtaining a preliminary distribution map of susceptibility property.However, Formula (1) may be used to update the distribution map ofsusceptibility property and the error limiting information. Other waysof obtaining the preliminary distribution map of susceptibility propertyare also applicable and within the scope of the present disclosure.Also, terms and iteration forms in Formula (1) may be modified oromitted. New terms or iteration forms may be introduced to Formula (1).Such modifications, deletions and/or variations are within theprotection scope of the present disclosure.

In 610, the updating unit 520 may perform an iterative process forupdating the preliminary distribution map of susceptibility property ofthe object and the preliminary error limiting information. In someembodiments, the iterative process may include at least one iteration.In some embodiments, the iterative process may include a plurality ofiterations. In each of the plurality of iterations, the error limitinginformation in the current iteration may be updated based on thedistribution map of susceptibility property in the preceding iteration;the distribution map of susceptibility property in the current iterationmay be updated based on the error limiting information updated in thecurrent iteration. Merely by way of example, the preliminarydistribution map of susceptibility property and the preliminary errorlimiting information may be treated as data in the zeroth iteration. Theerror limiting information in the first iteration may be determined (orupdated) according to the distribution map of susceptibility property inthe zeroth iteration (i.e., the preliminary distribution map ofsusceptibility property); the distribution map of susceptibilityproperty in the first iteration may be determined (or updated) accordingto the error limiting information determined in the first iteration.Similarly, the error limiting information in the second iteration may bedetermined according to the distribution map of susceptibility propertydetermined in the first iteration; the distribution map ofsusceptibility property in the second iteration may be determinedaccording to the error limiting information updated in the seconditeration, and so on. More descriptions regarding the updating of theerror limiting information and the distribution map of susceptibilityproperty may be found in, e.g., operations 804, 806 in FIG. 8, and thedescriptions thereof.

In some embodiments, after updating the distribution map ofsusceptibility property and the error limiting information in eachiteration, the updating unit 520 may determine whether a terminationcondition is satisfied. For example, the termination condition may besatisfied when a certain number or count of iterations have beenperformed in the iterative process. As another example, the terminationcondition may be satisfied when the difference between the field mapcorresponding to the updated distribution map of susceptibility propertyand the preliminary field map (also referred to as a convergence error)is less than a threshold. As a further example, the terminationcondition may be satisfied when the change of the distribution maps ofsusceptibility property in two or more successive iterations is lessthan a threshold. In some embodiments, the threshold may be determinedaccording to error information of an ideal distribution map ofsusceptibility property of the object. If the convergence error is lessthan or equal to the threshold, it may be determined that the accuracyof the distribution map of susceptibility property of the object issatisfied, and the iterative process may terminate. More descriptionsregarding the determination of whether the termination condition issatisfied may be found in, e.g., operation 808 in FIG. 8, and thedescriptions thereof.

In 612, the processing unit 530 may designate the updated distributionmap of susceptibility property of the object as a target distributionmap of susceptibility property of the object. For example, in responseto a determination that the termination condition is satisfied, or theplurality of iterations are completed, the iterative process mayterminate and the distribution map of susceptibility property in thecurrent (i.e., last) iteration may be designated as the targetdistribution map of susceptibility property of the object.

FIG. 7 is a flowchart of an exemplary process for determining apreliminary field map and/or a preliminary a priori map according tosome embodiments of the present disclosure. The process 700 may beexecuted by the processing device 140. For example, the process 700 maybe implemented as a set of instructions (e.g., an application) stored inthe storage, e.g., ROM 230, RAM 240, the storage device 150, the storage390, a storage device external to and accessible by the imaging system100. The processing device 140, the CPU 220, the CPU 340, the processingmodule 440, and/or the processing unit 530 may execute the set ofinstructions, and when executing the instructions, it may be configuredto perform the process 700. In some embodiments, operations 604 and 606in FIG. 6 may be performed according to at least a portion of theprocess 700.

In 702, the acquisition unit 510 may receive an MR signal from anobject. In some embodiments, the object or a region of interest (ROI,also referred to as a target site) of the object may be excited by an RFsignal. When the object or an ROI of the object relaxes, a spin signalor spin echo may be generated (and/or emitted) by the object (or the ROIof the object) and/or received by the imaging system 100. The spinsignal may be phase encoded to form a magnetic resonance (MR) signal. Insome embodiments, the MR signal received may be a complex signalincluding a real part and an imaginary part.

In 704, the processing unit 530 may generate a phase diagram from the MRsignal. In some embodiments, the phase diagram may be obtained byperforming an inverse tangent function on the ratio between the realpart and the imaginary part.

In 706, the processing unit 530 may determine a preliminary field mapbased on the phase diagram. In some embodiments, the preliminary fieldmap b(r) may be determined based on the phase diagram, as illustrated byFormula (2):ø(r)=Ø₀(r)+γB ₀ b(r)TE,  (2)where ø(r) denotes the phase diagram, Ø₀(r) denotes the initial phasethat relates to the structure of the RF signal transmitting coil(s), γdenotes the gyromagnetic ratio that is constant for the same nuclei andfield strength (e.g., the gyromagnetic ratio for a hydrogen proton in a3 Tesla magnetic field is 127.8 MHz), B₀ denotes the strength of themain magnetic field, b(r) denotes the field map (e.g., the preliminaryfield map), and TE denotes the echo time of the spin echo.

In 708, the processing unit 530 may generate a magnitude diagram fromthe MR signal. In some embodiments, the magnitude diagram may beobtained by taking a square root of the sum of squares of the real partand squares of the imaginary part of the MR signal. For example, for anMR signal S=a+i*b, the amplitude of the MR signal M=(a²+b²)^(0.5).

In 710, the processing unit 530 may determine a feature map based on themagnitude diagram. In some embodiments, properties of the magnitudediagram may be extracted to generate the feature map. The properties mayinclude the gradient property along one or more reference directions(e.g., the X direction, the Y direction, and/or the Z direction) of themagnitude diagram. The feature map may also be referred to as a gradientmap. For example, the magnitude diagram may include a plurality ofpixels, and each pixel of the magnitude diagram may have a pixel value.The feature map may include a plurality of pixels having a one-to-onecorrespondence with the plurality of pixels of the magnitude diagram.For each of the plurality of pixels in the feature map, the value of thepixel may be determined as an evaluation value between a pixel value ofa first pixel and a pixel value of a neighboring pixel adjacent to thefirst pixel of the magnitude diagram along a certain referencedirection. As used herein, two pixels are referred to as neighboringpixels if the two pixels are right next to each other and not spacedapart by another pixel along a certain direction.

In 712, the processing unit 530 may determine a preliminary a priori mapbased on the feature map. In some embodiments, the feature map obtainedin the operation 710 may be directly designated as the preliminary apriori map. Alternatively, in order to emphasize the boundaryinformation in the feature map, an image binarization may be performedon the feature map to generate the preliminary a priori map. Forinstance, the image binarization may include resetting pixel values ofpixels to zero or one by comparing the pixel values with a thresholdaccording to a rule. In some embodiments, if the pixel value of a pixelis lower than a threshold, the pixel value of the pixel may be reset tozero; if the pixel value of a pixel is higher than the threshold, thepixel value may be reset to one. If the pixel value of a pixel is equalto the threshold, the pixel value of the pixel may be reset to zero orone according to the rule. The rule may be set by the imaging system 100according to a default setting thereof, or provided by a user, etc. Insome embodiments, the threshold may be a relatively low number (e.g.,slightly lower than the gradient value of the pixels along innerboundaries) so that both an outer boundary (e.g., a boundary between theobject and air) and an inner boundary (e.g., a boundary betweendifferent types of tissue inside the object) may be retained.Alternatively, the threshold may be a relatively high number (e.g.,higher than the gradient value of the pixels along an inner boundary butlower than the gradient value of the pixels along an outer boundary) sothat only the outer boundary is retained.

In some embodiments, 708, 710 and 712 may be performed or omitted. If708, 710 and 712 are performed, the processing result may correspond tothe panel 1302; if 708, 710 and 712 are omitted, the processing resultmay correspond to the panel 1304.

FIG. 8 is a flowchart of an exemplary process for iteratively updatingthe distribution map of susceptibility property and/or error limitinginformation according to some embodiments of the present disclosure. Theprocess 800 may be executed by the processing device 140. For example,the process 800 may be implemented as a set of instructions (e.g., anapplication) stored in the storage, e.g., ROM 230, RAM 240, the storagedevice 150, the storage 390, a storage device external to and accessibleby the imaging system 100. The processing device 140, the CPU 220, theCPU 340, the processing module 440, and/or the processing unit 530 mayexecute the set of instructions, and when executing the instructions, itmay be configured to perform the process 800. In some embodiments,operation 610 in FIG. 6 may be performed according to at least a portionof the process 800.

In 802, the acquisition unit 510 may obtain a preliminary distributionmap of susceptibility property of the object and preliminary errorlimiting information (e.g., the preliminary a priori information and/orthe preliminary artifact estimation information). In some embodiments,the preliminary distribution map of susceptibility property of theobject may be generated by an operation described in operation 608 andthe preliminary error limiting information may be generated by anoperation described in operation 712. Merely by way of example, thepreliminary distribution map of susceptibility property and thepreliminary error limiting information may be treated as data in thezeroth iteration.

In 804, the updating unit 520 may update error limiting information inthe current iteration based on a distribution map of susceptibilityproperty in the preceding iteration. For example, in the firstiteration, the updating unit 520 may update error limiting informationin the first iteration based on a distribution map of susceptibilityproperty in the zeroth iteration (i.e., the preliminary distribution mapof susceptibility property). Similarly, the error limiting informationin the second iteration may be determined according to the distributionmap of susceptibility property determined in the first iteration.

In a case that the error limiting information includes the a prioriinformation (or the a priori map), more descriptions regarding theupdating of the error limiting information may be found in, e.g.,operations 904 and 906, and operations 908 and 910 in FIG. 9 and thedescriptions thereof (showing two different ways of updating the apriori information based on the distribution map of susceptibilityproperty). In a case that the error limiting information include theartifact estimation information (or an artifact estimation map, astreaking artifact estimation map), more descriptions regarding theupdating of the error limiting information may be found in, e.g.,operations 908, 912, and 914. In a case that the error limitinginformation includes both the a priori information and the artifactestimation information, the updating unit 520 may update the a prioriinformation and the artifact estimation information separately in anyorder or simultaneously.

In 806, the updating unit 520 may update the distribution map ofsusceptibility property in the current iteration based on the updatederror limiting information determined in the current iteration. Forexample, the distribution map of susceptibility property in the firstiteration may be determined (or updated) according to the error limitinginformation determined in the first iteration. Similarly, thedistribution map of susceptibility property in the second iteration maybe determined according to the error limiting information updated in thesecond iteration, and so on. In some embodiments, the distribution mapof susceptibility property may be updated based on the error limitinginformation using Formula (1).

As already discussed, Formula (1) describes an optimization process,e.g., seeking for a value χ to make the value of Formula (1) minimum orless than a threshold. In each of the plurality of iterations, Formula(1) may be used to determine a distribution map of susceptibilityproperty χ based on an updated difference map Δϕ (related to artifactestimation information) and an updated property mask P (related to the apriori information). The updated distribution map of susceptibilityproperty χ may then be used to update the difference map Δϕ and theproperty mask P. Merely by way of example, Formula (1) may be solved byany kinds of numerical optimization algorithm including, e.g., aconjugate gradient, a preconditioned conjugate gradient, quadraticminimization, split-Bregman, or the like, or any combination thereof.

It should be noted that Formula (1) includes both the a prioriinformation and the artifact estimation information. In a case that theerror limiting information only includes the a priori information, thethird part β∥F⁻¹D′FΔϕ∥₂ ² in Formula (1) may be omitted and Formula (1)may be reduced to:argmin_(χ) ∥F ⁻¹ DFχ−ϕ∥ ₂ ² +α∥P∇χ∥ ₁.  (3)

In a case that the error limiting information only includes the artifactestimation information, the term α∥P∇χ∥₁ in Formula (1) may be omittedand Formula (1) may be reduced to:argmin_(χ) ∥F ⁻¹ DFχ−ϕ∥ ₂ ² +β∥F ⁻¹ D′FΔϕ∥ ₂ ².  (4)

More descriptions regarding the updating of distribution map ofsusceptibility property based on the artifact estimation information maybe found in, e.g., FIG. 10 and the descriptions thereof.

In 808, the updating unit 520 may determine whether a terminationcondition is satisfied. For example, the termination condition may besatisfied when a certain number of count of iterations have beenperformed in the iterative process. As another example, the terminationcondition may be satisfied when the difference between field map relatedto the updated distribution map of susceptibility property and thepreliminary field map (also referred to as a convergence error) is lessthan a threshold. As a further example, the termination condition may besatisfied when the change of distribution maps of susceptibilityproperty in two or more successive iterations is less than a threshold.More descriptions regarding the determination as to whether thedifference between field map related to the updated distribution map ofsusceptibility property and the preliminary field map is less than athreshold may be found in, e.g., FIG. 11 and the descriptions thereof.In response to a determination that the termination condition is notsatisfied, the process 800 may return to operation 804, and the updatingunit 520 may perform a new iteration of updating the error limitinginformation and the distribution map of susceptibility property in 804and 806, respectively. In response to a determination that thetermination condition is satisfied, the process 800 may proceed tooperation 810.

In 810, the processing unit 530 may output the updated distribution mapof susceptibility property of the object and/or the updated errorlimiting information in the current iteration. The outputteddistribution map of susceptibility property of the object may bedesignated as a target distribution map of susceptibility property ofthe object.

In some embodiments, the process 800 optimizes the susceptibilitydistribution map through the iterative process. In each of theiterations in the iterative process, the a priori information and/or theartifact estimation information are dynamically updated so that the apriori information and/or the artifact estimation information graduallyapproach the true value, hence gradually reducing the convergence error(e.g., the evaluation value between the field map corresponding to themagnetic susceptibility map and the preliminary field map), andgradually improving the accuracy of the distribution map ofsusceptibility property of the object.

FIG. 9 is a flowchart of an exemplary process for updating the errorlimiting information according to some embodiments of the presentdisclosure. The process 900 may be executed by the processing device140. For example, the process 900 may be implemented as a set ofinstructions (e.g., an application) stored in the storage, e.g., ROM230, RAM 240, the storage device 150, the storage 390, a storage deviceexternal to and accessible by the imaging system 100. The processingdevice 140, the CPU 220, the CPU 340, the processing module 440, and/orthe processing unit 530 may execute the set of instructions, and whenexecuting the instructions, it may be configured to perform the process900. In some embodiments, operation 804 in FIG. 8 may be performedaccording to at least a portion of the process 900.

In 902, the acquisition unit 510 may obtain a distribution map ofsusceptibility property. In some embodiments, the distribution map ofsusceptibility property may be a preliminary distribution map ofsusceptibility property of the object or an updated distribution map ofsusceptibility property determined in a preceding iteration. Thedistribution map of susceptibility property determined during iterationsmay be updated based on operation 806.

In 904, the processing unit 530 may generate a feature map based on thedistribution map of susceptibility property. For example, the featuremap may be determined based on the distribution map of susceptibilityproperty of the object updated in the preceding iteration. As anotherexample, the feature map may be determined based on the distribution mapof susceptibility property in the preceding iteration and at least oneof the magnitude diagram, the phase diagram, the susceptibility weightedimaging (SWI) map, or the R2* (T2*) map. In some embodiments, propertiesof the distribution map of susceptibility property may be extracted togenerate the feature map. The properties may include the gradientproperty along one or more reference directions (e.g., the X direction,the Y direction, and/or the Z direction, etc.) of the distribution mapof susceptibility property. The feature map may also be referred to as agradient map. The distribution map of susceptibility property of theobject (e.g., updated in the current iteration) may include a pluralityof pixels, and each pixel of the distribution map of susceptibilityproperty of the object (e.g., updated in the current iteration) may havea pixel value. The feature map may include a plurality of pixels havinga one-to-one correspondence with the plurality of pixels of thedistribution map of susceptibility property. For each of the pluralityof pixels in the feature map, a evaluation value between a pixel valueof a first pixel and a pixel value of a neighboring pixel adjacent tothe first pixel along a certain reference direction in a precedingiteration may be determined as a value of a pixel in the feature map ofthe current iteration that corresponds to the first pixel in thedistribution map of susceptibility property of the object.

In 906, the updating unit 520 may update a priori information based onthe feature map. The a priori information may provide features on thespatial distribution of susceptibility of the object, such as thecontinuity of distribution within the same tissue, a boundary betweendifferent tissues, etc. For example, the a priori information mayinclude spatial smoothness of the susceptibility distribution of theobject, a spatial gradient of the susceptibility distribution of theobject, edge information of the susceptibility distribution of theobject, a weighted property mask, or the like, or any combinationthereof, corresponding to the target site of the object. In someembodiments, the updating unit 520 may generate a weighted property maskbased on the feature map and the weighted property mask may bedesignated as the a priori information of the current iteration. Theproperty mask may be an image including a plurality of pixels whosepixel values are, e.g., zeros and ones. For instance, values of pixelsin a smooth region (e.g., a region inside or outside a type of tissue ofthe object that has a relatively low gradient variance) may be zero,while values of pixels near a boundary having a relatively high gradientvariance may be one. In some embodiments, the values of pixels may beany value between zero and one. For example, the above-mentioned pixelsinside or outside a type of tissue with a low gradient variance may beany number between zero and one (e.g., 0.01, 0.1, 0.3) and theabove-mentioned pixels near a boundary may be any number between zeroand one (e.g., 0.99, 0.9, 0.7).

The a priori information may be dynamically generated and/or updatedbased on the intermediate distribution map of susceptibility propertyobtained during each iteration for obtaining the target distribution mapof susceptibility property. The updated a priori information may beapplied to the next iteration. In this way, the a priori information maygradually approach the spatial features of the ideal distribution map ofsusceptibility property, and the iterative convergence rate of theoptimization algorithm may also be increased. Furthermore, the accuracyof the results may be improved.

In 908, the processing unit 530 may generate an a priori field map basedon the distribution map of susceptibility property. In some embodiments,the a priori field map may be obtained by performing a forward magneticmodel function (e.g., using a forward dipole kernel in the forwardmagnetic model function) on the distribution map of susceptibilityproperty in a preceding iteration. In some embodiments, performing aforward magnetic model function on the distribution map ofsusceptibility property may also be referred to as performing a forwarddipole kernel function on the distribution map of susceptibilityproperty. In some embodiments, the distribution map of susceptibilityproperty in the current iteration may be obtained by performing aninverse magnetic model function (e.g., using an inverse dipole kernel inthe inverse magnetic model function) on the a priori field map in thecurrent iteration. In some embodiments, performing an inverse magneticmodel function on the a priori field map may also be referred to asperforming an inverse dipole kernel function on the a priori field map.The “a priori field map” may also be referred to as the “estimated fieldmap” as in the description of operations 912 and 914. It should be notedthat the a priori field map and the estimated field map are twodifferent names to refer to the same thing and they are usedinterchangeably in the present disclosure.

In some embodiments, the forward magnetic model function may beexpressed as:B(r)=B ₀ ·FT ⁻¹ {g(k)·χ(k)}.  (5)

The forward dipole kernel g(k) may be expressed as:

$\begin{matrix}{{g(k)} = \left\{ {\begin{matrix}{{\frac{1}{3} - \frac{k_{z}^{2}}{k_{x}^{2} + k_{y}^{2} + k_{z}^{2}}},} & {k \neq 0} \\{0,} & {k = 0}\end{matrix}.} \right.} & (6)\end{matrix}$

The inverse magnetic model function may be expressed as:(r)=FT ⁻¹ {g ⁻¹(k)·FT[ΔB(r)]}/B ₀,  (7)and the inverse dipole kernel g⁻¹(k) may be expressed as:

$\begin{matrix}{{g^{- 1}(k)} = \left\{ {\begin{matrix}{{3 \cdot \frac{k_{x}^{2} + k_{y}^{2} + k_{z}^{2}}{k_{x}^{2} + k_{y}^{2} - {2k_{z}^{2}}}},} & {k \neq 0} \\{0,} & {k = 0}\end{matrix},} \right.} & (8)\end{matrix}$where ΔB(r) denotes the field map; Δχ(r) denotes the distribution map ofsusceptibility property; FT and FT⁻¹ denotes Fourier transform matricesand inverse Fourier transform matrices, respectively; χ denotes themagnetic susceptibility distribution; k denotes a Fourier domaincoordinate vector; k_(x) denotes the coordinate in the X direction,k_(y) denotes the coordinate in the Y direction, k_(z) denotes thecoordinate in the Z direction; r denotes the coordinate vector in thespace domain; B₀ denotes the external magnetic field; and g and g⁻¹denote the forward magnetic dipole kernel transformation matrix and theinverse magnetic dipole kernel transformation matrix, respectively.

In 910, the updating unit 520 may update the a priori information basedon the a priori field map. In some embodiments, the a priori field mapmay be used as the a priori information in the current iteration.

In 912, the processing unit 530 may generate a difference field mapbased on the estimated field map and a preliminary field map. Merely byway of example, the pixel value of each pixel in the difference fieldmap may be calculated as the evaluation value between the pixel value ofa pixel of the preliminary field map and the pixel value of thecorresponding pixel of the estimated field map.

In 914, the updating unit 520 may update artifact estimation informationbased on the difference field map. In some embodiments, the artifactestimation information may be determined based on the difference fieldmap. For example, the forward dipole kernel D or a partial forwarddipole kernel D′ may be obtained. The artifact estimation informationmay be determined based on the forward dipole kernel D (or the partialforward dipole kernel D′) and the difference field map as illustrated inthe third part of Formula (1). The remnant radiation streaking artifactsmay be dynamically estimated during the iterative process of thedistribution map of susceptibility property, and may be uses as aregularization term in the iterative process. The operation may reducethe intensity of remnant radiation streaking artifacts to a negligiblelevel without affecting the results of other regularization terms,thereby reducing the artifacts and ensuring the accuracy of thecalculation results.

It should be noted that the above description of process 900 is merelyprovided for the purposes of illustration, and not intended to limit thescope of the present disclosure. For persons having ordinary skills inthe art, multiple variations and modifications may be made under theteachings of the present disclosure. However, those variations andmodifications do not depart from the scope of the present disclosure.For example, if the error limiting information only includes the apriori information, operations 912 and 914 may be omitted. As anotherexample, if the error limiting information only includes the artifactestimation information, operations 904, 906, and 910 may be omitted. Asa further example, operations 904 and 906 provide a first technique forupdating the a priori information, operations 908 and 910 provide asecond technique for updating the a priori information, and theoperations 904 and 906 and the operations 908 and 910 may be performedalternatively.

FIG. 10 is a flowchart of an exemplary process for updating thedistribution map of susceptibility property according to someembodiments of the present disclosure. The process 1000 may be executedby the processing device 140. For example, the process 1000 may beimplemented as a set of instructions (e.g., an application) stored inthe storage, e.g., ROM 230, RAM 240, the storage device 150, the storage390, a storage device external to and accessible by the imaging system100. The processing device 140, the CPU 220, the CPU 340, the processingmodule 440, and/or the processing unit 530 may execute the set ofinstructions, and when executing the instructions, it may be configuredto perform the process 1000. In some embodiments, if the error limitinginformation includes the artifact estimation information, operation 806in FIG. 8 may be performed according to at least a portion of theprocess 1000.

In 1002, the acquisition unit 510 may obtain updated artifact estimationinformation. For example, the artifact estimation information may beupdated based on the distribution map of susceptibility property and thepreliminary field map in operations 908, 912 and 914.

In 1004, the processing unit 530 may regularize the updated artifactestimation information. The regularization of the updated artifactestimation information may include introducing additional information inorder to solve an ill-posed problem (e.g., the singularity values in theinverse dipole kernel). For example, the artifact estimation information∥F⁻¹D′FΔϕ∥₂ ² in Formula (1) is regularized by a scaling factor β thatadjusts the weighting of the remnant streaking artifact regularizationterm.

In 1006, the updating unit 520 may update a distribution map ofsusceptibility property based on the regularized artifact estimationinformation. For example, the distribution map of susceptibilityproperty may be updated based on the regularized artifact estimationinformation β∥F⁻¹D′FΔϕ∥₂ ² and Formula (1) (or Formula (4)).

In 1008, the updating unit 520 may update the distribution map ofsusceptibility property based on the distribution map of susceptibilityproperty determined in the preceding iteration and the updated artifactestimation information determined in the current iteration. For example,the distribution map of susceptibility property in the current iterationmay be generated by subtracting the updated artifact estimationinformation determined in the current iteration from the distributionmap of susceptibility property in the preceding iteration.

It should be noted that the above description of process 1000 is merelyprovided for the purposes of illustration, and not intended to limit thescope of the present disclosure. For persons having ordinary skills inthe art, multiple variations and modifications may be made under theteachings of the present disclosure. However, those variations andmodifications do not depart from the scope of the present disclosure.For example, operations 1004 and 1006 provide a first technique forupdating the distribution map of susceptibility property based on theupdated artifact estimation information, operation 1008 provide a secondtechnique for updating the distribution map of susceptibility propertybased on the updated artifact estimation information, and the operations1004 and 1006 and the operation 1008 may be performed alternatively.

FIG. 11 is a flowchart of an exemplary process for determining whether atermination condition relating to the iteration process is satisfiedaccording to some embodiments of the present disclosure. The process1100 may be executed by the processing device 140. For example, theprocess 1100 may be implemented as a set of instructions (e.g., anapplication) stored in the storage, e.g., ROM 230, RAM 240, the storagedevice 150, the storage 390, a storage device external to and accessibleby the imaging system 100. The processing device 140, the CPU 220, theCPU 340, the processing module 440, and/or the processing unit 530 mayexecute the set of instructions, and when executing the instructions, itmay be configured to perform the process 1100. In some embodiments,operation 808 in FIG. 8 may be performed according to the process 1100.

In 1102, the acquisition unit 510 may obtain an updated distribution mapof susceptibility property. The distribution map of susceptibilityproperty may be updated as described in operation 806.

In 1104, the processing unit 530 may generate an estimated field map(also referred to as an a priori field map) based on the distributionmap of susceptibility property. In some embodiments, the estimated fieldmap may be generated in the current iteration by performing a partialforward dipole kernel function on the updated distribution map ofsusceptibility property of the object in the current iteration. Similarto operation 908, a forward dipole kernel function may be performed onthe distribution map of susceptibility property to generate theestimated field map.

In 1106, the processing unit 530 may generate a difference field mapbased on the estimated field map and a preliminary field map. Similar tooperation 912, the value of each pixel in the difference field map maybe determined based on a difference between the value of a pixel in theestimated field map and the value of the corresponding pixel in thepreliminary field map. As used herein, two pixels in two different mapsor images (e.g., a pixel in an estimated field map and a pixel in apreliminary field map) are considered corresponding to each other if thetwo pixels correspond to a same physical point in an object or in theimaging system 100.

In 1108, the processing unit 530 may determine an evaluation value (alsoreferred to as a convergence error) based on the difference field map.The evaluation value may represent the condition or the evaluationassociated with the difference field map. For example, the evaluationvalue may be the value of the L1 norm, L2 norm (referred in 608) or atotal variance. In some embodiments, the value of L2 norm may bedesignated as the evaluation value. The evaluation value may beexpressed as the square root of the sum of squares of values of all thepixels in the difference field map. For example, the evaluation valuemay be expressed as:

$\begin{matrix}{{{{Evaluation}\mspace{14mu}{value}} = \sqrt{\sum\limits_{{i = {1:X}};{j = {1:Y}};{k = {1:Z}}}\;\left( {{values}\mspace{14mu}{of}\mspace{14mu}{pixels}\mspace{14mu}\left( {x_{i},y_{j},z_{k}} \right)^{2}} \right.}},} & (9)\end{matrix}$wherein X denotes the size of the difference field map along the X axis;Y denotes the size of the difference field map along the Y axis; and Zdenotes the size of the difference field map along the Z axis.

In 1110, the processing unit 530 may determine whether the evaluationvalue is less than a threshold. In response to a determination that theevaluation value is less than the threshold, the process 1100 mayproceed to operation 1112; in response to a determination that theevaluation value exceeds the threshold, the process 1100 may proceed tooperation 1114. If the evaluation value is equal to the threshold, theprocess 1100 may proceed to operation 1112 or operation 1114 accordingto a rule. The rule may be set by the imaging system 100 according to adefault setting thereof, or provided by a user, etc.

In 1112, the processing unit 530 may determine that the terminationcondition is satisfied, i.e., the updating of the distribution map ofsusceptibility property may be completed. In some embodiments, inresponse to a determination that the termination condition is satisfied,the process 1100 may proceed to operation 810.

In 1114, the processing unit 530 may determine that the terminationcondition is not satisfied, i.e., the updating of the distribution mapof susceptibility property is not completed. In some embodiments, inresponse to a determination that the termination condition is notsatisfied, the process 1100 may return to operation 804.

FIG. 12 is a schematic diagram illustrating an exemplary idealdistribution map of susceptibility property, an exemplary field map, andan exemplary reconstructed distribution map of susceptibility propertywith artifacts according to some embodiments of the present disclosure.In some embodiments, the image 1202 may be an unknown but a targetdistribution map of susceptibility property of a cylinder. The image1204 may be an acquired field map of the cylinder in image 1202. Theimage 1206 may be a distribution map of susceptibility property that isdirectly reconstructed from the image 1204. It should be noted that theimage 1206 contains apparent artifacts. The goal of the techniquesdisclosed in the present application is to generate the image 1202 basedon the image 1204.

EXAMPLES

The following examples are shown for illustration purposes and notintended to limit the scope of the present disclosure.

FIG. 13 is a schematic diagram illustrating exemplary schemes ofiteratively updating a priori information according to some embodimentsof the present disclosure. As shown in FIG. 13, an example of a specificexperimental result of iteratively (or dynamically) updating a prioriinformation is shown. The a priori information in this example is aweighted property mask obtained from the spatial gradient features(e.g., the gradient variance) of the image. The image spatial gradientwas calculated independently in the three directions of the X direction,the Y direction, and the Z direction as illustrated in FIG. 1, and threecorresponding gradient maps Gx, Gy, and Gz were obtained. The gradientmaps Gx, Gy, and Gz were extracted as a mask (Gx), a mask (Gy), and amask (Gz), respectively. In the present example, pixels were setaccording to a thresholding technique. When the value of a certain pixelin the gradient map(s) exceeded the threshold, the pixel value of pixelsin the corresponding weighted property mask was set to 0 (white);otherwise, the pixel value of pixels in the corresponding weightedproperty mask was set to 1 (black).

Panel 1302 shows mask Gx, mask Gy, and mask Gz determined in the firstiteration through the third iteration, and the ninth iteration of theiterative process. The initial weighted property mask (i.e., before thefirst iteration) was determined based on the magnitude diagram. In eachiteration of the iterative process of updating the distribution map ofsusceptibility property, the weighted property mask was updated based onthe distribution map of susceptibility property in the precedingiteration and used for updating the distribution map of susceptibilityproperty in the current iteration. In some embodiments, the panel 1302may correspond to the processing result when the steps 708, 710 and 712are performed.

Panel 1304 shows mask Gx, mask Gy, and mask Gz determined in the firstiteration through the third iteration, and the ninth iteration of theiterative process. The initial weighted property mask was determinedbased on a brain mask, e.g., the initial weighted property mask does notcontain useful a priori knowledge. In each iteration of the iterativeprocess of updating the distribution map of susceptibility property, theweighted property mask was updated based on the distribution map ofsusceptibility property in the preceding iteration and used for updatingthe distribution map of susceptibility property in the currentiteration. In some embodiments, the panel 1304 may correspond to theprocessing result when the steps 708, 710 and 712 are not performed.

Panel 1306 shows mask Gx, mask Gy, and mask Gz determined in the firstiteration through the third iteration, and the ninth iteration of theiterative process. The initial weighted property mask was determinedbased on the magnitude diagram. In each iteration of the iterativeprocess of updating the distribution map of susceptibility property, theweighted property mask was updated based on both the distribution map ofsusceptibility property in the preceding iteration and the initialweighted property mask (e.g., by a superposition operation). Hence, thedifference between panel 1302 and panel 1306 is that the weightedproperty mask in panel 1306 contains both the spatial a prioriinformation of the intermediate susceptibility map and the spatial apriori information of the magnitude diagram, while the panel 1302 onlycontains the spatial a priori information of the intermediatesusceptibility map.

It should be noted that by combining both the intermediate distributionmap of susceptibility property and the initial weighted property mask todynamically (or iteratively) update the a priori information as shown inpanel 1306, a clearer property mask may be obtained. Therefore, aclearer and more accurate distribution map of susceptibility propertymay be generated based on the clearer property mask. It should be notedthat all the techniques disclosed in generating masks illustrated inpanels 1302, 1304, and 1306 may be implemented by the imaging system 100and are within the protection scope of the present disclosure.

FIG. 14 is a schematic diagram illustrating exemplary computer simulateddemonstrations of artifact of various amount in the distribution map ofsusceptibility property according to some embodiments of the presentdisclosure. The left column 1402 is an ideal distribution map ofsusceptibility property generated by computer simulation. The middlecolumn 1404 is a difference map between the ideal distribution map ofsusceptibility property in the left column 1402 and the targetdistribution map of susceptibility property generated according to theprocesses disclosed in, e.g., FIGS. 6-11 of the present disclosure(e.g., dynamically updating the a priori information, the remnantstreaking artifact estimation information and the distribution map ofsusceptibility property). The right column 1406 is a difference mapbetween the ideal distribution map of susceptibility property in theleft column 1402 and the distribution map of susceptibility propertygenerated by the conventional QSM technique. The gray levels of thepixels of the ideal distribution map of susceptibility property in theleft column 1402 may be in a range of −0.3˜0.3 part per million (ppm),in which the highest gray level may be 0.3 ppm or even higher, while thelowest gray level may be −0.3 ppm or even lower. Each of the pixels ofthe difference maps in the middle column 1404 and the right column 1406is presented in difference in brightness indicating its value. Forexample, the darkest pixel in the difference maps may correspond to avalue of −0.05 part per million (ppm), and the brightest pixel in thedifference maps may correspond to a value of 0.05 ppm. It can be seenfrom FIG. 14 that the values of pixels of the difference map in themiddle column 1404 are much closer to zero (bright pixels are lessbright and dark pixels are less dark) than the difference map in theright column 1406. This may indicate that the distribution map ofsusceptibility property determined by the present technique described inthe present disclosure is more accurate with less convergence error thanthe conventional QSM technique.

FIG. 15 is a schematic diagram illustrating exemplary human brain imagesdemonstrating artifacts of various amount in the distribution maps ofsusceptibility property according to some embodiments of the presentdisclosure. Column 1502 illustrates three distribution maps ofsusceptibility property generated according to the technique disclosedin, e.g., FIGS. 6-11 of the present disclosure (e.g., dynamicallyupdating the a priori information, the remnant streaking artifactestimation information and the distribution map of susceptibilityproperty). Column 1504 illustrates three distribution maps ofsusceptibility property generated according to the conventional QSMtechnique. Column 1506 illustrates three difference maps between column1502 and column 1504. Each row of columns 1502-1506 illustrates athree-dimensional human brain image obtained based on the same scandata. As shown in FIG. 15, the images in column 1502 have a relativelylow level of streaking artifacts, indicating that the techniquedescribed in the present disclosure is effective in eliminating orsuppressing artifacts compared with the conventional QSM technique.

Having thus described the basic concepts, it may be rather apparent tothose skilled in the art after reading this detailed disclosure that theforegoing detailed disclosure is intended to be presented by way ofexample only and is not limiting. Various alterations, improvements, andmodifications may occur and are intended to those skilled in the art,though not expressly stated herein. These alterations, improvements, andmodifications are intended to be suggested by this disclosure, and arewithin the spirit and scope of the exemplary embodiments of thisdisclosure.

Moreover, certain terminology has been used to describe embodiments ofthe present disclosure. For example, the terms “one embodiment,” “anembodiment,” and/or “some embodiments” mean that a specific feature,structure or characteristic described in connection with the embodimentis included in at least one embodiment of the present disclosure.Therefore, it is emphasized and should be appreciated that two or morereferences to “an embodiment” or “one embodiment” or “an alternativeembodiment” in various portions of this specification are notnecessarily all referring to the same embodiment. Furthermore, thespecific features, structures or characteristics may be combined assuitable in one or more embodiments of the present disclosure.

Further, it will be appreciated by one skilled in the art, aspects ofthe present disclosure may be illustrated and described herein in any ofa number of patentable classes or context including any new and usefulprocess, machine, manufacture, or composition of matter, or any new anduseful improvement thereof. Accordingly, aspects of the presentdisclosure may be implemented entirely hardware, entirely software(including firmware, resident software, micro-code, etc.) or combiningsoftware and hardware implementation that may all generally be referredto herein as a “unit,” “module,” or “system.” Furthermore, aspects ofthe present disclosure may take the form of a computer program productembodied in one or more computer readable media having computer readableprogram code embodied thereon.

A computer readable signal medium may include a propagated data signalwith computer readable program code embodied therein, for example, inbaseband or as part of a carrier wave. Such a propagated signal may takeany of a variety of forms, including electro-magnetic, optical, or thelike, or any suitable combination thereof. A computer readable signalmedium may be any computer readable medium that is not a computerreadable storage medium and that may communicate, propagate, ortransport a program for use by or in connection with an instructionexecution system, apparatus, or device. Program code embodied on acomputer readable signal medium may be transmitted using any appropriatemedium, including wireless, wireline, optical fiber cable, RF, or thelike, or any suitable combination of the foregoing.

Computer program code for carrying out operations for aspects of thepresent disclosure may be written in any combination of one or moreprogramming languages, including an object oriented programming languagesuch as Java, Scala, Smalltalk, Eiffel, JADE, Emerald, C++, C#, VB. NET,Python or the like, conventional procedural programming languages, suchas the “C” programming language, Visual Basic, Fortran 2003, Perl, COBOL2002, PHP, ABAP, dynamic programming languages such as Python, Ruby andGroovy, or other programming languages. The program code may executeentirely on the user's computer, partly on the user's computer, as astand-alone software package, partly on the user's computer and partlyon a remote computer or entirely on the remote computer or server. Inthe latter scenario, the remote computer may be connected to the user'scomputer through any type of network, including a local area network(LAN) or a wide area network (WAN), or the connection may be made to anexternal computer (for example, through the Internet using an InternetService Provider) or in a cloud computing environment or offered as aservice such as a Software as a Service (SaaS).

Furthermore, the recited order of processing elements or sequences, orthe use of numbers, letters, or other designations therefore, is notintended to limit the claimed processes and methods to any order exceptas may be specified in the claims. Although the above disclosurediscusses through various examples what is currently considered to be avariety of useful embodiments of the disclosure, it is to be understoodthat such detail is solely for that purpose, and that the appendedclaims are not limited to the disclosed embodiments, but, on thecontrary, are intended to cover modifications and equivalentarrangements that are within the spirit and scope of the disclosedembodiments. For example, although the implementation of variouscomponents described above may be embodied in a hardware device, it mayalso be implemented as a software only solution, e.g., an installationon an existing server or mobile device.

Similarly, it should be appreciated that in the foregoing description ofembodiments of the present disclosure, various features are sometimesgrouped together in a single embodiment, figure, or description thereoffor the purpose of streamlining the disclosure aiding in theunderstanding of one or more of the various embodiments. This method ofdisclosure, however, is not to be interpreted as reflecting an intentionthat the claimed subject matter requires more features than areexpressly recited in each claim. Rather, claimed subject matter may liein less than all features of a single foregoing disclosed embodiment.

What is claimed is:
 1. A method implemented on a computing device havingat least one storage device storing a set of instructions fordetermining a distribution map of susceptibility property of an objectand at least one processor in communication with the at least onestorage device, the method comprising: obtaining a phase diagramcorresponding to a magnetic resonance (MR) signal of the object;determining a preliminary field map based on the phase diagram;obtaining preliminary error limiting information associated with thepreliminary field map; determining, based on the preliminary field mapand the preliminary error limiting information, a preliminarydistribution map of susceptibility property of the object; performing aniteration process including at least one iteration, wherein each of theat least one iteration comprises: updating error limiting information ina current iteration based on a distribution map of susceptibilityproperty of the object updated in a preceding iteration; updating thedistribution map of susceptibility property of the object in the currentiteration based on the error limiting information updated in the currentiteration; and performing a comparison based on the updated distributionmap in the current iteration and the preliminary field map; anddesignating the updated distribution map of susceptibility property ofthe object as a target distribution map of susceptibility property ofthe object.
 2. The method of claim 1, wherein the performing acomparison based on the updated distribution map in the currentiteration and the preliminary field map comprises: generating anestimated field map based on the updated distribution map ofsusceptibility property of the object in the current iteration;comparing the estimated field map with the preliminary field map; anddetermining whether the updated distribution map in the currentiteration satisfies a condition based on a result of the comparisonbetween the estimated field map and the preliminary field map.
 3. Themethod of claim 2, wherein the comparing the estimated field map withthe preliminary field map comprises: determining a difference field mapbetween the estimated field map and the preliminary field map;determining an evaluation value based on the difference field map; andthe determining whether the updated distribution map in the currentiteration satisfies a condition comprising: comparing the evaluationvalue with a threshold.
 4. The method of claim 1, wherein the errorlimiting information comprises an a priori map, and the obtainingpreliminary error limiting information comprises: determining at leastone of a magnitude diagram, a phase diagram, a susceptibility weightedimaging (SWI) map, or a R2* (T2*) map based on the MR signal; anddetermining a preliminary a priori map based on the at least one of themagnitude diagram, the phase diagram, the susceptibility weightedimaging (SWI) map, or the R2* (T2*) map.
 5. The method of claim 1,wherein the error limiting information comprises an a priori map, andthe updating error limiting information in a current iteration based onthe distribution map of susceptibility property of the object updated ina preceding iteration comprises: determining a feature map in thecurrent iteration based at least in part on the distribution map ofsusceptibility property of the object updated in the precedingiteration; and updating the a priori map based on the feature map. 6.The method of claim 5, wherein the distribution map of susceptibilityproperty of the object updated in the current iteration includes aplurality of pixels, each pixel of the distribution map ofsusceptibility property of the object in the current iteration having apixel value, the feature map includes a plurality of pixels having aone-to-one correspondence with the plurality of pixels of thedistribution map of susceptibility property of the object, and thedetermining a feature map in the current iteration based on thedistribution map of susceptibility property of the object updated in thepreceding iteration comprises: for each of the plurality of pixels inthe feature map, determining an evaluation value between a pixel valueof a first pixel and a pixel value of a second pixel adjacent to thefirst pixel in the distribution map of susceptibility property of theobject in the preceding iteration; and designating the evaluation valueas a pixel value of a pixel in the feature map of the current iterationthat corresponds to the first pixel in the distribution map ofsusceptibility property of the object.
 7. The method of claim 6, whereinthe determining an evaluation value between a pixel value of a firstpixel and a pixel value of a second pixel next to the first pixel in thedistribution map of susceptibility property of the object in thepreceding iteration comprises: determining the evaluation value betweenthe pixel value of the first pixel and the pixel value of the secondpixel next to the first pixel in a reference direction, the referencedirection including at least one of an X direction, a Y direction, or aZ direction.
 8. The method of claim 1, wherein the error limitinginformation comprises an a priori map, and the updating error limitinginformation in a current iteration based on the distribution map ofsusceptibility property of the object in a preceding iterationcomprises: determining an a priori field map in the current iterationbased on the distribution map of susceptibility property of the objectin the preceding iteration; and updating the a priori map in the currentiteration based on the a priori field map in the current iteration. 9.The method of claim 8, wherein determining an a priori field map in thecurrent iteration based on the distribution map of susceptibilityproperty of the object in the preceding iteration comprises: determiningthe a priori field map in the current iteration by performing a forwarddipole kernel function on the distribution map of susceptibilityproperty of the object in the preceding iteration.
 10. The method ofclaim 1, wherein the error limiting information comprises an artifactestimation map, and the updating error limiting information in a currentiteration based on the distribution map of susceptibility property ofthe object in a preceding iteration comprises: generating an estimatedfield map in the current iteration based on the updated distribution mapof susceptibility property of the object in the current iteration;determining a difference field map between the estimated field map inthe current iteration and the preliminary field map; and updating theartifact estimation map in the current iteration based on the differencefield map.
 11. The method of claim 10, wherein generating an estimatedfield map in the current iteration based on the updated distribution mapof susceptibility property of the object in the current iterationcomprises: generating the estimated field map in the current iterationby performing a partial forward dipole kernel function on the updateddistribution map of susceptibility property of the object in the currentiteration.
 12. The method of claim 10, wherein the updating thedistribution map of susceptibility property of the object in the currentiteration based on the error limiting information updated in the currentiteration comprises: generating the updated distribution map ofsusceptibility property of the object in the current iteration bysubtracting an updated artifact estimation map in the current iterationfrom the distribution map of susceptibility property of the object inthe preceding iteration; or generating the updated distribution map ofsusceptibility property of the object in the current iteration by usinga regularization term relating to an updated artifact estimation map inthe iteration process.
 13. The method of claim 1, wherein the errorlimiting information includes an a priori map and an artifact estimationmap, and the updating the distribution map of susceptibility property ofthe object in the current iteration based on the error limitinginformation updated in the current iteration comprises: generating theupdated distribution map of susceptibility property of the object in thecurrent iteration based on the a priori map in the current iteration,the artifact estimation map in the current iteration and the preliminaryfield map.
 14. The method of claim 1, wherein the iteration process isperformed based on at least one of a conjugate gradient, apreconditioned conjugate gradient, quadratic minimization, orsplit-Bregman.
 15. The method of claim 1, wherein the error limitinginformation includes an a priori map relating to a spatial smoothness, aspatial gradient, edge information, or a weighted property mask of thedistribution map of susceptibility property of the object.
 16. A systemcomprising: at least one storage including a set of instructions orprograms for determining a distribution map of susceptibility propertyof an object; at least one processor configured to communicate with theat least one storage, wherein when executing the set of instructions orprograms, the at least one processor is directed to perform operationsincluding: obtaining a phase diagram corresponding to a magneticresonance (MR) signal of the object; determining a preliminary field mapbased on the phase diagram; obtaining preliminary error limitinginformation associated with the preliminary field map; determining,based on the preliminary field map and the preliminary error limitinginformation, a preliminary distribution map of susceptibility propertyof the object; performing an iteration process including at least oneiteration, wherein in each of the at least one iteration includes:updating error limiting information in a current iteration based on adistribution map of susceptibility property of the object updated in apreceding iteration; updating the distribution map of susceptibilityproperty of the object in the current iteration based on the errorlimiting information updated in the current iteration; and performing acomparison based on the updated distribution map in the currentiteration and the preliminary field map; and designating the updateddistribution map of susceptibility property of the object as a targetdistribution map of susceptibility property of the object.
 17. Thesystem of claim 16, wherein to perform a comparison based on the updateddistribution map in the current iteration and the preliminary field map,the at least one processor is directed to perform additional operationsincluding: generating an estimated field map based on the updateddistribution map of susceptibility property of the object in the currentiteration; comparing the estimated field map with the preliminary fieldmap; and determining whether the updated distribution map in the currentiteration satisfies a condition based on a result of the comparisonbetween the estimated field map and the preliminary field map.
 18. Thesystem of claim 16, wherein the error limiting information comprises ana priori map, and to obtain preliminary error limiting information, theat least one processor is directed to perform additional operationsincluding: determining at least one of a magnitude diagram, a phasediagram, a susceptibility weighted imaging (SWI) map, or a R2* (T2*) mapbased on the MR signal; and determining a preliminary a priori map basedon the at least one of the magnitude diagram, the phase diagram, thesusceptibility weighted imaging (SWI) map, or the R2* (T2*) map.
 19. Thesystem of claim 16, wherein the error limiting information comprises anartifact estimation map, and to update error limiting information in acurrent iteration based on the distribution map of susceptibilityproperty of the object in a preceding iteration, the at least oneprocessor is directed to perform additional operations including:generating an estimated field map in the current iteration based on theupdated distribution map of susceptibility property of the object in thecurrent iteration; determining a difference field map between theestimated field map in the current iteration and the preliminary fieldmap; and updating the artifact estimation map in the current iterationbased on the difference field map.
 20. A non-transitory computerreadable medium comprising executable instructions or programs that,when executed by at least one processor, cause the at least oneprocessor to effectuate a method comprising: obtaining a phase diagramcorresponding to a magnetic resonance (MR) signal of the object;determining a preliminary field map based on the phase diagram;obtaining preliminary error limiting information associated with thepreliminary field map; determining, based on the preliminary field mapand the preliminary error limiting information, a preliminarydistribution map of susceptibility property of the object; performing aniteration process including at least one iteration, wherein each of theat least one iteration includes: updating error limiting information ina current iteration based on a distribution map of susceptibilityproperty of the object updated in a preceding iteration; updating thedistribution map of susceptibility property of the object in the currentiteration based on the error limiting information updated in the currentiteration; and performing a comparison based on the updated distributionmap in the current iteration and the preliminary field map; anddesignating the updated distribution map of susceptibility property ofthe object as a target distribution map of susceptibility property ofthe object.